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Predicting  MoM  Error  Currents  by  Inverse 
Application  of  Residual  E-Fields 

Andre  P.C.  Fourie,  Derek  C.  Nitch,  and  Alan  R.  Clark 


Abstract —  This  paper  presents  a  methodology  to  predict 
a  posteriori  the  error  associated  with  a  Method-of-Moments 
solution.  The  discussion  is  limited  to  a  one  dimensional 
pulse  basis  function  wire-based  implementation,  but  is  eas¬ 
ily  extended.  A  Formulation  for  Error  Prediction  based  on 
the  relationship  between  the  error  in  the  boundary  condi¬ 
tions  and  the  error  in  the  solution  is  presented,  and  vali¬ 
dated  by  an  over- segmented  problem.  The  formulation  is 
then  used  in  a  normally-segmented  solution  to  predict  the 
error  by  means  of  a  linear  interpolation  of  the  calculated  cur¬ 
rent  which  results  in  a  smoother  boundary  condition  error. 
The  results  show  that  this  normally-segmented  methodology 
predicts  the  error  current  within  5%  of  the  “accurate”  error 
current  obtained  by  a  20:1  oversegmentation  of  the  problem. 
Further  work  needs  to  be  performed  to  extend  this  to  the 
multidimensional  case,  although  no  technical  difficulties  are 
expected  with  this. 

I.  Introduction 

In  a  method  of  moments  (MoM)  solution  to  an  electro¬ 
magnetic  problem,  the  structure  currents  are  calculated  to 
ensure  that  the  boundary  conditions  are  satisfied  in  some 
sense,  usually  at  specific  match  points  (point  matching)  or 
over  specific  domains  in  an  average,  or  weighted  average, 
manner. 

The  boundary  condition  to  be  satisfied — when  using  a 
MoM  for  perfectly  conducting  wires — is  that  the  total  tan¬ 
gential  E-field  should  be  zero  at  all  points  on  the  wire.  A 
current  which  ensures  zero  tangential  E-fields  in  a  contin¬ 
uous  sense  on  all  wires  would  be  accurate  in  accordance 
with  the  uniqueness  theorem. 

Many  researchers  have  recognized  that  the  accuracy  of 
a  method  must  be  related  to  how  well  the  boundary  con¬ 
ditions  are  met  by  a  specific  solution  (Hsaio  &  Kleinman 
1996,  Meyer  &  Davidson  1996).  The  exact  relationship 
between  the  error  in  the  boundary  conditions  and  the  er¬ 
ror  in  the  observable  quantities  has  not  been  defined  and 
has  also  not  been  used  to  gain  an  estimate  of  solution  (or 
method)  accuracy.  This  study  investigates  the  possibility 
of  getting  such  a  relationship  and  then  formulates  approx¬ 
imations  which  allow  its  incorporation  in  MoM  (and  other 
methods)  with  the  lowest  computational  overhead. 

II.  Formulation  for  Error  Prediction  (FEP) 

The  obvious  relationship  between  the  error  in  the  bound¬ 
ary  conditions  and  the  error  in  the  solution  is  quite 
clear  in  the  following  formulation  for  perfectly  conducting 
wires  (Thiele  1973): 

Suppose  we  have  obtained  a  current  solution,  /(s),  where 
I(s)=Ia(s)+Ie(s)  (1) 

The  authors  are  with  the  Department  of  Electrical  Engineering, 
University  of  the  Witwatersrand,  Johannesburg,  Wits,  2050  South 
Africa. 


and  Ia(s)  is  the  accurate  solution  and  Ie(s)  is  the  error 
associated  with  the  solution  I(s). 

If  we  have  the  continuous  (and  assumed  accurate  for  the 
moment)  relationship  between  the  current  and  the  excita¬ 
tion  then 

Lopla(s)  =  V(s)  (2) 

where  Lop  is  the  linear  operator  (Pocklington’s  equation  is 
an  example  for  wire  problems)  which  defines  the  excitation 
tangential  to  the  wire,  V  (s).  Pocklington’s  equation  is  usu¬ 
ally  solved  using  the  Method-of-Moments.  In  the  case  of  a 
point  matching  solution,  this  results  in  a  tangential  E-field 
error  of  zero  only  at  the  match  points.  At  any  other  point 
on  the  structure  there  will  be  a  tangential  E-field  error. 

If  pulse  weighting  functions  were  used  in  the  MoM  solu¬ 
tion,  then  it  is  likely  that  there  will  be  a  tangential  E-field 
error  at  any  point  on  the  structure — the  tangential  E-field 
error  will  be  zero  only  in  the  average  sense. 

The  tangential  E-field  errors  found  on  the  structure  are 
defined  to  be  the  residuals.  The  tangential  E-field  from  the 
inaccurate  current,  I(s),  is: 

LopI(s)  =  V(s)  +  Ve(s)  (3) 

where  Ve(s )  is  the  tangential  E-field  residual  (or  error)  on 
the  wires.  Combining  equations  (1)  and  (3),  and  from  the 
linearity  of  the  operator  and  superposition  it  follows  that 

Lople(s)  =  ve(s)  (4) 

The  inverse  of  equation  (4)  will  allow  the  error  current  Ie  (s) 
to  be  obtained  from  the  E-field  residual  (error),  Ve(s). 

Ie(s)  =  L~pVe(s)  (5) 

We  shall  call  the  above  mathematical  development  the  For¬ 
mulation  for  Error  Prediction  (FEP)  for  later  reference. 
Once  the  error  currents  have  been  found,  the  errors  in 
secondary  parameters  such  as  radiation  pattern  or  input 
impedance  are  easily  quantified.  This  may  be  done,  for  ex¬ 
ample,  by  using  the  error  current  to  obtain  near/far  fields 
at  the  same  points  as  the  full  solution  values  were  com¬ 
puted.  This  will  produce  “error  fields”  which  places  a 
bound  on  the  error  in  the  computed  field  values.  Errors 
in  other  parameters  such  as  input  impedance,  etc.  may  be 
obtained  in  a  similar  manner. 

This  somewhat  trivial  proof  indicates  that  the  error  in 
any  solution  may,  in  principle,  be  calculated  exactly  pro¬ 
vided  that: 

•  The  residual  E-field  over  the  structure  can  be  calculated 
continuously  and  accurately. 
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•  An  accurate  continuous  operator  Lop  and  its  inverse  ex¬ 
ist. 

In  practice  the  following  obstacles  to  using  this  result  are 
apparent: 

•  Calculation  of  the  residual  E-field  over  the  structure  can 
be  computationally  time  intensive. 

•  The  operator  is  normally  discretized  as  an  interaction 
matrix  and  hence  is  neither  continuous  nor  accurate. 

III.  Using  the  FEP  to  obtain  accurate  error 

CURRENTS. 

The  graphs  which  follow  illustrate  an  investigation  using 
the  FEP  for  a  specific  problem.  A  dipole  antenna  of  length 
l  =  0.5 A  and  with  radius,  a  =  0.005 A  was  excited  with  an 
incident  E-field  of  1  V/m. 

A  simple  MoM  program  was  written  which  uses  pulse 
basis  functions  and  point  matching.  The  convergence  of  the 
method  is  shown  in  Figure  1.  The  method  clearly  converges 
to  a  stable  answer  when  around  100  to  200  segments  are 
used  for  the  dipole.  The  Formulation  for  Error  Prediction 


Fig.  1.  The  feed  current  magnitude  for  a  variation  in  number  of 
segments  of  a  dipole  with  l  =  0.5A  and  a  =  0.005A  using  pulse 
basis  functions  and  point  matching. 

(FEP)  is  illustrated  by  considering  the  current  from  the 
200  segment  problem  to  be  the  “accurate”  current.  Two 
“inaccurate”  currents  were  obtained  by  solving  the  same 
problem  with  only  10  and  20  segments. 

The  residual  E-fields  from  these  two  inaccurate  solu¬ 
tions  were  then  obtained  at  200  points  across  the  dipole 
(i.e.  at  many  more  points  than  the  match  points).  Figure 
2  shows  the  residual  E-field  magnitude  for  a  10  segment 
solution  and  a  20  segment  solution,  which  clearly  demon¬ 
strates  that  the  boundary  conditions  are  exactly  met  at  the 
match  points,  but  also  shows  a  significant  error  (50 V/m  for 
a  lV/m  excitation)  at  other  points.  The  residual  E-fields 
in  Figure  2  can  be  applied  as  an  excitation  vector  on  the 
same  geometry  divided  into  200  segments  to  yield  the  error 
current  in  accordance  with  the  FEP  as  stated  in  equations 
(4)  and  (5).  Figure  3  shows  this  graphically.  It  is  clear 


Fig.  2.  The  magnitude  of  the  residual  E-fields  obtained  for  a  10  and 
20  segment  dipole  with  l  —  0.5A  and  a  =  0.005A  using  pulse  basis 
functions  and  point  matching. 


from  these  figures  that  sum  of  the  inaccurate  current  (ob¬ 
tained  from  the  10  segment  solution)  and  the  error  current 
(obtained  by  applying  the  residual  E-fields  via  FEP)  is 
practically  equal  to  the  “accurate”  current  obtained  from 
the  200  segment  case.  The  200  segment  inverse  interaction 
matrix  hence  represents  an  accurate  and  pseudo-continuous 
inverse  operator,  Lop,  and  the  200  point  discretization  of 
the  residual  E-field  a  pseudo-continuous  representation  on 
the  error,  Ee(s).  Figure  3  illustrates  the  principle  behind 


Fig.  3.  The  current  magnitude  to  prove  the  FEP  for  a  10  segment 
dipole  with  t  =  0.5A  and  a  =  0.005A  using  pulse  basis  functions 
and  point  matching. 


error  prediction  by  inverse  application  of  the  residual  E- 
fields  as  excitations  to  produce  an  error  current.  The  com¬ 
putational  effort  does  not  render  it  suitable  as  a  general 
tool  for  error  prediction — oversegmenting  a  problem  by  a 
factor  of  20  in  order  to  obtain  an  error  estimate  is  clearly 
unproductive!  The  technique  may  however,  be  useful  as 
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Fig.  4.  The  error  current  associated  with  NEC2  solutions  for  a  2,  10 
segment  dipole  ,  as  compared  to  the  10  segment  MoM  solution 
using  pulses  with  i  =  0.5A  and  a  =  0.005A. 


a  research  tool;  different  basis  and/or  weighting  functions, 
for  instance,  may  be  compared  in  terms  of  absolute  error 
currents.  Such  an  exercise  was  performed  to  compare  the 
NEC2  (Burke  &  Poggio  1981)  basis  functions  with  pulse 
basis  functions.  The  same  dipole  antenna  was  simulated 
using  NEC2  with  2  and  10  segments  and  the  error  currents 
obtained  using  the  same  procedure,  by  applying  the  tan¬ 
gential  E-field  errors  as  an  excitation  to  an  oversegmented 
wire  of  the  same  dimensions.  The  error  currents  associated 
with  these  two  cases  are  shown  in  Figure  4,  together  with 
the  10  segment  MoM  problem,  presented  earlier.  Com¬ 
paring  these  currents  indicates  that  NEC2  basis  functions 
achieve  roughly  the  same  accuracy  for  2  segments  as  pulse 
basis  functions  do  for  10  segments.  This  may  be  expected 
since  NEC2  uses  a  superpositioned  sine,  cosine  and  a  con¬ 
stant  term  as  its  basis  function,  which  is  quite  suitable  for 
approximating  the  currents  on  a  dipole. 

IV.  Using  FEP  without  oversegmentation 

The  previous  section  showed  that  the  FEP  may  be  used 
to  obtain  accurate  error  currents  with  major  computa¬ 
tional  overheads  associated  with  overdiscretizing  the  prob¬ 
lem.  We  now  attempt  to  use  the  same  FEP  without  the 
requirement  to  increase  the  problem  discretization.  The 
main  problem  with  the  FEP  is  that  although  it  is  possible 
to  calculate  the  residual  accurately,  only  N  excitations  may 
be  used  to  excite  the  error  currents.  It  is  hence  necessary 
to  reduce  the  residual  over  each  segment  to  a  single  value. 
The  first  approximation  to  the  accurate  FEP  was  to  obtain 
the  residual  E-fields  at  100  points  on  the  dipole  and  use  the 
average  value  of  ten  samples  applied  as  Ee(s )  on  the  origi¬ 
nal  10  segment  problem1.  Results  from  using  this  method 

1Note  that  using  an  “average  value”  is  meaningful  for  a  point 
matching  MoM  (using  Dirac  delta  weighting  functions),  however,  if 
pulse  weighting  functions  are  used,  the  error  residual  computed  would 
be  zero.  In  the  case  of  the  pulse  weighting  function,  one  sample  per 
segment  may  result  in  a  more  appropriate  segment  error.  It  stands 


were  not  satisfactory,  since  the  large  oscillations  in  residual 
E-field  evident  in  Figure  2  require  a  finer  discretization  in 
order  to  obtain  a  reasonable  estimate — which  defeats  the 
objective.  It  should  be  noted  that  the  estimate  was  not  in¬ 
accurate  as  a  result  of  insufficient  field  point  samples  when 
averaging,  since  increasing  the  samples  to  200  and  400  did 
not  improve  matters. 


Fig.  5.  The  current  obtained  for  a  10  segment  dipole  with  l  =  0.5 A 
and  a  =  0.005 A  using  pulse  basis  functions  and  point  matching 
and  a  linear  interpolation  to  the  original  pulse  current. 


Fig.  6.  The  tangential  E-field  magnitude  from  a  10  segment  dipole 
with  t  =  0.5 A  and  a  =  0.005 A  using  pulse  basis  functions  and 
point  matching  as  well  as  the  tangential  E-field  from  a  linear 
approximation  to  the  original  pulse  current 

We  recognized  that  the  large  oscillations  in  the  resid¬ 
ual  E-fields  were  mainly  due  to  the  pulse  basis  functions 
which  are  discontinuous  at  segment  boundaries.  Rather 
than  solving  the  problem  with  different  basis  functions  to 

to  reason  that  the  method  used  to  obtain  the  segment  residual  must 
employ  a  function  other  than  that  of  the  actual  testing  functions  used 
in  the  MoM  solution. 
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Fig.  7.  Error  current  from  applying  the  averaged  E-tangential  fields 
obtained  from  the  pulse  currents  and  the  linearly  interpolated 
currents  as  an  excitation  to  the  10  segment  problem.  The  accu¬ 
rate  error  current  obtained  earlier  is  also  shown  for  comparison. 

ensure  smooth  residual  E-field  behaviour,  we  performed 
curve  fitting  to  the  actual  pulse  current  solution  before  cal¬ 
culating  the  residual  E-fields  via  FEP.  Figure  5  shows  the 
original  pulse  currents  obtained  from  the  MoM  together 
with  a  linear  approximation  to  these  currents,  which  clearly 
offer  a  smoother  current  solution.  Figure  6  compares  the 
residual  E-field  due  to  the  original  pulse  currents  as  well 
as  to  the  linearly  interpolated  currents  (using  only  10  seg¬ 
ments).  The  residuals  from  the  interpolated  currents  are 
clearly  much  less  oscillatory  than  those  resulting  from  the 
actual  pulse  solution. 

The  average  (over  a  segment)  of  the  residual  E-fields  as¬ 
sociated  with  the  linearly  interpolated  current  shown  in 
Figure  6  was  applied  to  the  original  10  segment  problem 
to  obtain  error  currents  without  increasing  problem  dis¬ 
cretization.  Figure  7  shows  the  effectiveness  of  this  ap¬ 
proach.  Using  the  average  residual  E-fields  obtained  from 
the  linearly  interpolated  currents  provides  a  good  estimate 
of  the  accuracy  of  the  solution. 

It  should  be  noted  that,  in  general,  the  linear  interpola¬ 
tion  of  the  pulse  basis  function  solution  will  not  necessarily 
yield  a  more  accurate  current  representation.  (It  was  also 
not  the  purpose.)  The  interpolated  current  was  merely 
used  to  produce  better  behaved  (less  oscillatory)  residuals. 
These  residuals  can  be  averaged  over  one  segment  to  yield 
a  residual  E-field  vector  of  the  same  order  as  the  problem 
discretization  (10  values  in  our  example).  The  FEP  can, 
hence,  be  used  with  the  same  size  matrix,  and  if  matrix  fac¬ 
torization  is  used,  without  much  additional  computational 
effort. 

V.  Conclusion 

The  FEP  is  definitely  suitable  for  research  purposes  when 
applying  it  “accurately”  by  increasing  the  problem  dis¬ 
cretization.  An  error  prediction  scheme  for  estimating 
the  errors  associated  with  run-of-the-mill  simulations  us¬ 


ing  MoM,  however,  would  not  be  useful  if  problem  order  is 
increased.  This  is  mainly  due  to  memory  limitations  (since 
matrix  storage  space  is  proportional  to  N 2)  and  to  com¬ 
putational  limitations  (since  execution  time  is  proportional 
to  iV3).  The  latter  part  of  this  limited  study  shows  that 
a  reasonable  estimate  may  be  obtained  without  using  finer 
problem  discretization  provided  that  the  behaviour  of  the 
E-field  residuals  are  smooth  over  a  segment.  This  is  true 
even  if  a  discontinuous  current  due  to  simple  basis  func¬ 
tions  is  artificially  smoothed  by  interpolation  to  achieve 
better  behaviour  of  residual  E-fields. 

It  seems  to  us  that  the  relationship  between  residual  E- 
field  smoothness  and  required  problem  discretization  must 
obey  normal  sampling  criteria  for  the  FEP  to  provide  rea¬ 
sonable  error  estimates.  Problem  discretization  should  al¬ 
low  the  residual  E-field  curve  shape  to  be  retained  by  the 
error  excitation  vector  for  good  estimates.  Residual  E- 
field  variations  along  a  wire  with  a  period  of  the  order  of 
the  segment  length  (0.05 A  for  10  segments)  are  the  result 
of  non-physical  behaviour  of  currents  at  segment  bound¬ 
aries;  it  should  hence,  in  principle,  always  be  possible  to 
“smooth”  current  behaviour  to  obtain  slower  variations  in 
residuals  (changes  in  the  order  of  a  wavelength  should  be 
more  realistic). 

The  interpolation  approach  also  illustrates  a  very  inter¬ 
esting  aspect  of  the  behaviour  of  the  residual  E-fields:  the 
large,  fast,  oscillations  are  purely  associated  with  local, 
non-physical  current  behaviour.  Smoothing  the  current 
results  in  much  lower  variations  in  the  residual  E-fields, 
which  render  them  more  suitable  for  error  prediction.  The 
interpolated  currents  are  not,  however,  inherently  more 
accurate — at  least  in  terms  of  errors  in  overall  magnitude — 
but  they  do  seem  to  be  a  more  natural  representation  of 
current  shape,  which  one  would  expect.  In  the  cases  in¬ 
vestigated  here,  current  interpolation  was  done  in  order 
to  obtain  error  estimates  with  lower  computational  over¬ 
heads.  The  opposite  is  clearly  also  possible:  currents  may 
be  obtained  using  some  crude  basis  functions  and  can  then 
be  smoothed  a  posteriori  (or  altered)  while  using  the  FEP 
method  to  measure  whether  a  more  accurate  answer  was 
obtained. 

Naturally,  errors  in  the  secondary  parameters  such  as 
input  impedance  or  radiation  pattern,  etc.  can  be  derived 
from  the  error  currents  using  the  existing  relationships  be¬ 
tween  currents  and  these  parameters — allowing  for  “error 
bars”  to  be  placed  on  these  parameters. 
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Abstract  -  In  an  effort  to  mathematically  validate  the 
convergence  properties  of  various  surface  current-iterative 
methods,  the  magnetic  field  integral  equation  is  analyzed  for 
its  contraction  mapping  properties.  The  analysis  is 
performed  first  on  the  general  integral  operators  and  then  on 
the  matrices  representing  the  discrete  forms  of  the  integral 
operators  associated  with  the  different  iterative  methods . 
The  contraction  mapping  properties  are  determined  by 
investigating  the  spectral  radius  of  each  linear  operator. 
Conditions  for  the  verification  and  validation  of  these 
iterative  methods  are  provided ,  along  with  mathematical 
checks  for  the  existence  of  spurious  modes  and  the  existence 
of  internal  resonance. 

1.  Introduction 

Throughout  the  past  decade,  a  variety  of  surface  current- 
iterative  methods  [1-6]  have  been  proposed  to  solve 
electromagnetic  scattering  and  radiation  problems  involving  the 
Magnetic  Field  Integral  Equation  (MFIE)  and  related  integral 
equations.  Most  of  these  methods  were  developed  as 
alternatives  to  the  traditional  Method  of  Moments  (MoM)  [7,8] 
solution,  which  involve  the  inversion  of  a  dense  {most  often 
complex)  matrix.  While  each  of  the  methods  showed  some 
computational  advantages,  each  also  introduced  a  new  issue  of 
concern:  that  of  convergence.  For  the  most  part,  the 
convergence  of  these  iterative  methods  has  been  shown  by 
numerical  example,  but  little  mathematical  analysis  has  been 
offered  for  the  understanding  of  why  these  methods  converge  to 
the  correct  solution.  Oftentimes  convergence  is  simply  claimed 
when  the  delta  change  in  iterations  falls  below  some  pre-set 
limit.  In  the  following  paper,  the  MFIE  is  examined  to  explain 
why  these  surface  current-iterative  methods  have  been 
successfully  iterated  to  convergence. 

In  previous  work  Kaye  et  al[  1]  found  that  iterative  solutions  to 
the  MFIE  for  the  surface  current  would  always  yield 
convergence  for  Perfectly  Electrically  Conducting  (PEC) 
bodies  which  where  divided  into  two  parts:  a  side  illuminated 
by  the  incident  field,  and  a  shadowed  side.  Coupled  integral 


equations  were  written  for  each  side.  The  iterative  process 
started  on  the  illuminated  side  by  iterating  from  the  Physical 
Optics  (PO)  [9]  current  to  an  improved  value.  The  process 
then  went  to  the  shadowed  side  where  the  improved  current  on 
the  illuminated  side  was  used  with  an  initial  estimate  {most 
often  zero)  for  the  shadowed  side  current  to  obtain  an  improved 
shadowed  side  current.  The  process  continued  by  returning  to 
the  illuminated  side  to  obtain  further  improvement  for  the 
illuminated  side  current  and  then  returned  to  the  shadowed  side 
to  obtain  further  improvement  to  the  shadowed  current.  The 
process  was  halted  when  the  delta  change  in  the  solution  fell 
below  some  pre-set  limit. 

Reuster  &  Thiele  [2]  found  that  the  same  iterative  procedure 
could  be  used  for  PEC  cavities  where  the  aperture  of  the  cavity 
was  treated  as  the  illuminated  side  and  the  cavity  walls  were 
treated  as  the  shadowed  side.  The  iterative  process  began  by 
taking  the  aperture  field  to  be  that  produced  by  an  external 
plane  wave  illumination.  The  aperture  field  was  then  used  to 
find  the  total  magnetic  field,  via  iteration,  along  the  cavity 
walls.  The  field  along  the  cavity  walls  was  then  used  to  update 
the  total  field  at  the  aperture  via  iteration.  The  improved 
knowledge  of  the  aperture  field  was  then  used  to  obtain  a 
further  update  of  the  field  along  the  cavity  walls,  which  in  turn 
was  used  to  update  the  aperture  field.  Again,  the  process  was 
halted  when  the  delta  change  in  the  solution  fell  below  some 
pre-set  limit. 

Obelleiro-Basteiro  et  al  [3]  demonstrated  a  method  similar  to 
the  iterative  method  proposed  in  [2]  where  the  PO 
approximation  was  used  to  simplify  the  iterative  procedure. 
Numerical  results  are  presented  which  demonstrate  the 
convergence  and  accuracy  of  the  method,  but  no  mathematical 
analysis  is  provided  as  to  why  the  method  converges.  Collins 
&  Skinner  [4]  demonstrated  an  iterative  method  for  calculating 
the  scattering  from  perturbed  circular  dielectric  cylinders  by 
using  equivalent  currents  along  the  perimeter  of  the  cylinder. 
Their  paper  alludes  to  the  mathematical  properties  of  the 
method's  convergence,  but  no  mathematical  analysis  is 
presented.  Reuster,  et  al.  [5]  showed  by  numerical  example 
that  convergence  could  be  obtained  by  directly  iterating  the 


1054-4887  ©  1999  ACES 


REUSTER,  THIELE,  ELOE:  ITERATION  OF  SURFACE  CURRENTS  AND  THE  MAGNETIC  I.E. 


77 


entire  currents  on  a  PEC  body  if  the  PO  approximation  was 
used  as  an  initial  condition;  however,  little  mathematical 
analysis  was  provided.  Finally,  Hodges  &  Rahmat-Samii  [6] 
present  an  advanced  iterative  method  for  large  PEC  bodies 
consisting  of  both  wires  and  closed  surfaces.  Their  iterative 
method  involves  the  Electric  Field  Integral  Equation  (EFIE)  as 
well  as  the  MFTE  and  the  PO  approximation.  Again,  no 
mathematical  analysis  is  provided  as  to  why  the  method 
converges,  and  convergence  is  determined  when  the  delta 
change  in  the  solution  falls  below  some  pre-set  limit. 

Essentially,  all  of  these  surface  current-iterative  methods  are 
fixed-point  iterative  problems  [10-12]  where  the  MFIE  or  a 
related  integral  equation  is  manipulated  ( usually  by  observing 
physical  characteristics  associated  with  the  particular 
problem )  to  develop  an  iterative  scheme  which  is  “hopefully”  a 
contraction  mapping.  The  establishment  of  a  contraction 
mapping  is  the  unique  feature  that  guarantees  that  the  iterative 
method  will  converge  in  a  monotonic  mean-square  sense.  In 
the  work  that  follows,  two  formulations  of  the  magnetic  field 
integral  equation  (Maue's  formulation  and  the  Total  Field 
formulation )  are  analyzed  for  their  potential  contraction 
mapping  characteristics.  The  effect  of  the  discrete  form 
operator  size  on  the  contraction  mapping  properties  of  each 
formulation  is  studied,  and  a  mathematical  check  is  presented 
for  the  existence  of  spurious  modes  and  the  existence  of 
internal  resonance.  Finally,  a  general  iterative  scheme,  which 
involves  subdividing  the  integral  operator  for  the 
creation/insurance  of  contraction  mappings,  is  presented.  This 
general  iterative  scheme  is  related  to  the  surface  current- 
iterative  methods  [1-6]  that  were  discussed  earlier. 


2.  Contraction  Mapping  Analysis  of  the  MFIE 

For  simplicity  a  2D  Transverse  Electric  (TE)  PEC  scattering 
problem,  as  shown  in  Figure  1,  is  chosen  as  a  baseline  model 
for  analysis.  This  particular  scattering  problem  allows  the 
normally  complex  3-dimensional  vector  integral  equation  to  be 
reduced  to  a  2-dimensional  scalar  integral  equation  which  still 
maintains  all  of  the  properties  associated  with  the  general 
vector  integral  equation.  Hence,  no  loss  in  generality  is 
experienced  in  working  with  the  reduced  integral  equation. 
The  derivation  of  both  Maue's  equation  and  the  Total  Field 
equation  proceeds  as  follows. 


-1™-  - y\rfITz( r' )hf>( P\r-r/\)cos<pdl  =  0 
r  —>r  4  *  ,  and 

/2  H z  (  r  ^  is  the  principle  value  of  the  integral.  From  (1),  both 
Maue’s  formulation  and  the  Total  Field  formulation  for  the 
MFIE  may  be  obtained  as  follows. 


Figure  1  -  2D  TE  PEC  Scattering  Problem 


Maue’s  Formulation.  Maue's  formulation  is  obtained  from  (1) 
by  subtracting  the  integral  operator  from  both  sides  of  (1)  and 
then  scaling  (1)  by  a  factor  of  two. 

Htz( r)=-2H{( P \r-r'\)cosQdl  ^ 


Total  Field  Formulation.  The  Total  Field  formulation  is 
obtained  from  (1)  by  subtracting  the  integral  operator  from  both 
sides  of  (1)  and  then  adding  the  principle  value  to  both  sides  of 
(1). 

HTz(r)=-Hlz(r)+ Jc H Tz ( r ' ) h(?H P  \  r  -  r ' || cos Qdl 


+  jHTz(r) 


(3) 


Note  that  the  only  significant  difference  between  Maue's 
formulation  and  the  Total  Field  formulation  is  the  location  of 
the  principle  value  of  the  integral.  It  will  be  shown  that  the 
location  of  the  principle  value  has  a  major  effect  on  the 
contraction  mapping  properties  of  the  MFIE. 


For  2D  electromagnetic  radiation  and  scattering  problems, 
where  the  magnetic  field  is  parallel  to  the  geometry  of  interest, 
the  MFIE  may  be  stated  as  [1,13]: 

-H[(~r)=  ^\cHTz(r')hf(P\r-r'\)cos<l>dl  +  ^HTz(r) 

Where  C  is  the  simple  closed  curve  (in  the  x-y  plane) 
describing  the  PEC  body  of  interest. 


Integral  equations  (2)  and  (3)  may  be  solved  using  fixed-point 
iteration,  provided  that  these  integral  formulations  are 
contraction  mappings.  By  definition,  contraction  mappings  may 
be  iterated  directly  to  a  unique  solution  with  a  guarantee  of 
mean  squared-monotonic  convergence  [10-12].  This 
convergence  is  independent  of  the  initial  guess  used  to  initiate 
the  iterative  process.  The  general  test  for  a  contraction 
mapping  proceeds  as  follows. 
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Let  ¥ i  and  ¥2  represent  two  approximations  for  ¥  in  an 
equation  of  the  following  form,  where  L  is  any  linear  operator 
and  F  is  a  constant  forcing  function 

y/  =  L(i//J+F.  (4) 


HTz(rm)=-2nUrm) 

HTz(rn )h(i2)( P\rm'r„  \)COS(Pmn 


(9) 


Then  (4)  is  a  contraction  mapping  and  may  be  solved  directly 
using  fixed-point  iteration  (with  the  guarantee  of  mean  squared 
monotonic  convergence)  if  and  only  if 

\UW2  ~Vi)\-k\V2  ~¥i\ *  0<k<1  (5) 

Applying  the  above  test  to  Maue's  formulation  yields  (6). 

\HTz(r)2-HTz(rh\> 

4  Jr br (~r'h-HTz( r'h\h(P( P\  r-r  | jcos <j)dl 
2  ^  (6) 


After  simplification,  using  the  triangle  inequality  theorem,  the 
contraction-mapping  test  can  be  obtained  for  Maue's 
formulation  (7). 


Jc  |  hp(  P I  r  -  r  | )  cos  0 1  dl  <  2- 


(7) 


Similarly,  the  contraction-mapping  test  can  be  obtained  for  the 
Total  Field  formulation  (8). 


jc\hp(  P\r-r'\)cos<p\dl<j 


(8) 


Note  that  (7)  and  (8)  are  purely  functions  of  the  problem's 
geometry  and  the  wavelength  of  interest  ( P  =  ).  The 

requirements  for  convergence  are  independent  of  the  initial 
conditions  used  in  the  iterative  process.  Also  note  that  the 
convergence  condition  for  Maue’s  equation  is  more  difficult  to 
satisfy  than  the  convergence  condition  for  the  Total  Field 
equation.  This  is  a  direct  result  of  the  removal  of  the  principle 
value  from  the  integral  operator  in  Maue's  formulation.  While 
both  (7)  and  (8)  may  be  shown  true  for  a  particular  geometry, 
there  does  not  exist  a  general  proof  for  arbitrary  PEC  scattering 
bodies.  Hence,  further  investigation  into  the  contraction 
mapping  properties  of  both  Maue's  formulation  and  the  Total 
Field  formulation  is  continued  on  a  discretized/numerical  level. 


3.  Contraction  Mapping  Analysis  of  the  Discretized  MFIE 

The  inability  to  work  with  (7)  and  (8)  results  in  the  need  to 
approximate  the  integral  operators  in  (2)  and  (3)  with  systems 
of  linear  equations.  Applying  a  pulse-basis  point-patching 
MoM  expansion  [7]  to  (2)  and  (3)  results  in  the  following 
approximate  expansions  for  Maue's  formulation  and  the  Total 
Field  formulation.  Similar  analysis  can  be  performed  for  more 
advanced  MoM  expansions  [8]. 


HTzGm)=-H[(rm) 

4  2  (IQ) 


Expressing  (9)  and  (10)  in  matrix  from  yields  (11)  and  (12), 


hIGJ 
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Note  that  the  principle  value  of  the  integral  appears  as  the 
diagonal  elements  of  the  matrix  equation  associated  with  the 
Total  Field  formulation  (12)  while  the  principle  value  of  the 
integral  is  incorporated  into  the  forcing  function  of  the  matrix 
equation  associated  with  Maue's  formulation  (11).  This  is  a 
direct  result  of  the  removal  of  the  principle  value  from  the 
integral  operator  in  Maue's  formulation. 

It  is  now  possible  to  analyze  (11)  and  (12)  for  their  contraction 
mapping  properties  in  a  fashion  similar  to  the  analysis 
performed  in  Section  2.0.  For  matrix  equations  in  the  form  of 
(11)  and  (12),  the  matrix  equation  represents  a  contraction 
mapping  if  and  only  if  the  spectral  radius  of  the  matrix  operator 
P  is  less  than  one  [10,1 1].  The  spectral  radius  of  a  matrix  is 
defined  as  the  magnitude  of  the  largest  eigenvalue  of  the 
matrix.  In  general,  eigenvalues  are  the  most  difficult  of  all 
matrix  characteristics  to  compute,  and  for  this  reason  an 
advanced  matrix  analysis  package  was  utilized  to  perform  the 
following  study.  The  matrix  analysis  package  is  LAPACK,  and 
was  developed  for  Linear  Operator  analysis.  This  software 
package  is  in  the  public  domain,  and  may  be  obtained  via  the 
Internet. 

While  the  magnitudes  of  the  eigenvalues  associated  with  a 
given  MoM  expansion  are  typically  functions  of  the  number  of 
elements  used  in  the  expansion,  it  may  be  shown  that  the 
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eigenfunctions  of  a  given  MoM  expansion  converge  in  a 
fashion  similar  to  the  convergence  of  the  MoM  solution. 
Figures  2a  and  2b  show  the  magnitude  of  the  complex 
eigenvalues  for  a  0.5  A  radius  PEC  circular  cylinder  using  20, 
30,  and  40  basis  functions. 


Eigenvalue  Convergence  (Maue) 


radius  =  0.5  wavelengths 


Eigenvalue  Convergence  (Total) 


radius  =  0.5  wavelengths 


Eigenvalue  Number 


(b) 

Figure  2  -  Eigenvalue  Convergence 

For  both  Maue’s  formulation  and  the  Total  field  formulation, 
convergence  occurs  at  approximately  30  basis  functions  (which 
is  approximately  8  basis  functions  per  wavelength).  This  is 
consistent  with  MoM  expansions  for  a  smooth  surface  [7].  It 
should  be  noted  that  no  new  eigenvalues  occur  above  30  basis 
functions,  which  equates  to  no  new  information  being  gained 
by  increasing  the  number  of  basis  functions  being  used.  Also, 
note  that  the  eigenvalues  for  the  Total  field  equation  tend 
toward  0.5  (for  large  eigenvalue  numbers)  and  the  eigenvalues 
for  Maue’s  formulation  tend  toward  0.0  (for  large  eigenvalue 


numbers).  These  two  limiting  eigenvalues  directly  correspond 
to  the  diagonal  values  of  their  respective  integral  expansion 
(11)  and  (12). 

Figures  3a-3d  show  plots  of  the  complex  eigenvalue  versus  the 
eigenvalue  number  for  2-dimensional  circular  cylinders  with 
radii  of  0.60985  A  and  0.38275  A  ,  respectively.  These  radii 
were  chosen  because  they  correspond  to  the  TE01  and  TM01 
cutoff  frequencies  for  circular  waveguides  [13].  Circular 
cylinders  with  these  radii  are  historically  difficult  to  solve 
numerically  because  of  internal  resonance  problems.  It  should 
be  noted  that  in  both  cases  Maue's  equation  has  eigenvalues 
with  magnitudes  greater  than  one,  and  the  Total  Field  equation 
has  eigenvalues  with  magnitudes  less  than  or  equal  to  one.  This 
implies,  for  these  particular  cases,  direct  iteration  of  Maue's 
equation  will  not  yield  convergence  and  direct  iteration  of  the 
Total  Field  equation  may  yield  convergence  (the  largest 
eigenvalue  has  to  be  strictly  less  than  one  to  guarantee 
convergence).  In  addition,  for  radii  that  correspond  to  the 
cutoff  frequencies  of  the  TEOn  and  TMOn  circular  waveguide 
modes,  both  the  Total  Field  equation  and  Maue’s  equation  will 
have  eigenvalues  of  value  one.  It  is  these  eigenvalues  (of  value 
one)  which  are  responsible  for  the  historically  documented 
resonance  problems.  It  also  should  be  noted  that  Total  Field 
equation  has  a  zero  eigenvalue  for  the  0.60985  A  case  (figure 
3b).  This  zero  eigenvalue  can  cause  a  spurious  mode  to  appear 
in  the  solution.  While  eigenvalues  of  value  zero  are  not 
generally  problems  for  direct  iterative  solutions,  eigenvalues  of 
values  one  or  greater  are  problems,  and  lead  to  divergent 
solutions.  It  is  this  problem  of  eigenvalues  of  value  one  or 
greater  that  the  previous  surface  current-iterative  methods  [1-6] 
are  indirectly  solving.  Section  4  presents  a  general  method  for 
directly  dealing  with  the  problem  of  eigenvalues  of  value  one  or 
greater.  Furthermore,  it  may  be  shown  that  the  eigenvalues  for 
Maue’s  equation  and  those  for  the  Total  Field  equation  are 
simply  related. 

Maue’s  Formulation 

Circular  Cylinder  -  Radius  0.60985 


Complex  Eigenvalues  Complex  Eigenvalues  Complex  Eigenvalues 
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Total  Field  Formulation 

Circular  Cylinder  -  Radius  0.60985 


Maue’s  Formulation 


Circular  Cylinder  -  Radius  0.38275 


Total  Field  Formulation 

Circular  Cylinder  -  Radius  0.38275 


Figure  3  -  Complex  Eigenvalue  Analysis 


Write  Maue's  equation  from  (2)  in  operator  form  as 
HTz  =  Lm(H)-2H{ 

and  the  Total  Field  equation  from  (4)  as 
Htz  =  Lt(H)-h[ 

> 

where  the  operators  ll  and  Lm  are  related  by 

Lr  =  +  ~Lm 

and  [  I  ]  is  the  identity  matrix. 


(13) 

(14) 

(15) 


Let  ?.n,vm  be  an  eigenvalue  and  eigenvector  pair  for  Lm ,  then 
vm  is  an  eigenvector  for  LT  corresponding  to  an  eigenvalue 

^  ”  Yl  ^  +  ^  .  To  see  this,  assume 


LV)=1V, 


niv \ 


It  follows  that 
LT(vm)=^[l  +  Lmlvm 

=\[vm+Lm(vm)]=^[l+Xm]vm= XTvm 
2  2  .  (17) 

By  inspection  of  the  curves  in  Figures  2  and  3,  it  may  be  seen 

that  the  above  relationship  is  true  for  large  eigenvalue  numbers. 

By  noting  the  complex  nature  of  the  eigenvalues  in  Figure  3  it 

may  also  be  seen  that  the  above  relationship  is  true  for  the  small 

eigenvalue  numbers. 


4.  Methods  of  Insuring  a  Contraction  Mapping 

As  can  be  seen  from  Figures  2  and  3,  not  all  of  the  2D  TE  PEC 
scattering  problems  contain  a  spectral  radius  that  is  less  than 
one.  In  general,  the  spectral  radius  of  the  matrix  associated 
with  the  Total  Field  formulation  is  less  the  spectral  radius  of  the 
matrix  associated  with  Maue’s  formulation.  However, 
geometrical  cases  still  exist  where  the  spectral  radius  of  the 
matrix  associated  with  the  Total  Field  formulation  is  greater 
than  or  equal  to  one  (figure  3b  and  3d).  For  these  situations, 
the  resulting  system  of  linear  equations  can  not  be  solved  using 
a  direct  implementation  of  the  fixed-point  iterative  method,  and 
attempts  to  do  so  will  result  in  a  rapid  divergence  in  the 
solution  vector.  For  this  case,  where  the  spectral  radius  of  the 
matrix  associated  with  the  particular  magnetic  field  formulation 
is  greater  than  one,  a  modified  version  of  the  fixed-point 
iterative  method  must  be  applied  if  a  contraction  mapping  like- 
iterative  solution  is  to  be  obtained. 

In  their  work  [1,14-16]  Thiele  et  al.  showed  that  it  was  possible 
to  apply  the  method  of  fixed-point  iteration  to  a  2D  TE  PEC 
scattering  problem  by  subdividing  the  scattering  body  into  two 
pieces  which  they  defined  as  the  illuminated  side  and  the 
shadow  side.  While  their  work  was  restricted  to  a  single 
subdivision  of  Maue's  formulation,  their  method  is  shown  here 
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to  be  completely  general  and  may  be  applied  to  any  of  the 
magnetic  field  formulations  for  any  number  of  arbitrary 
subdivisions  provided  the  largest  eigenvalue  of  each 
subdivision  is  less  than  unity. 

We  can  express  the  hybrid  iterative  method  (HIM)  [1]  as  (18) 
where  the  sub-matrix  l  Lu  ]  is  associated  with  the  illuminated 
side,  the  sub-matrix  /  L22]  is  associated  with  the  shadowed 
side,  and  the  sub-matrices  IL12 1  and  [L21]  are  associated 
with  the  coupling  between  the  illuminated  and  shadowed  sides. 
The  iterative  method  used  in  [1]  to  solve  (18)  is  represent  in 
(19). 


L Lu )  [L12] 

Hl(rJ 

_  [L21]  \ L22 \ . 

htzGJ 

2  HlCrJ 

1 - 1 

0 

1 _ 1 

= 

■k,j  [or 
[0]  w. 

HlCrJ 

2H[CrJ 

[0] 

HlCrJ 

= 

■  [/]  [or 

.[zj  [/]. 

HlCrJ 

- 

[0] 

2H,zCrm) 

HlCrJ 

= 

■[/]  [0]- 
.[0]  u. 

HlCrJ 

- 

[0] 

2H[Crm) 

HlGj 

= 

m  tzj- 
jo]  [/]. 

HlCj 

2H{Crm) 

[0] 

- 1 

£' 

0 

= 

HlCrJ 

j 

_ 

. 

This  four  step  iterative  method  was  found  to  always  converge 
without  internal  resonance  problems,  although  a  formal  proof 
was  never  given  [1].  If  we  examine  the  eigenvalues  of  the 
resulting  sub-matrices,  it  is  apparent  why  the  method  converged 
and  why  it  had  no  internal  resonance  problems.  Figures  4a-d 
show  the  eigenvalues  for  the  HIM  matrices 
[  Liil  J  L12]  J  L21]  and  [  L22]  of  a  circular  cylinder  with 
radius  0.5  wavelengths.  The  eigenvalues  for  each  of  the  four 


matrices  are  all  less  than  one;  hence,  each  of  the  four  sub¬ 
iterations  is  itself  a  contraction  mapping.  Since  the  four  sub¬ 
iterations  are  nested  in  a  linear  fashion,  then  by  the  triangle 
inequality  theorem  the  total  iterative  procedure  is  itself  a 
contraction  mapping.  Thus,  the  system  of  linear  equations  may 
be  iterated  directly  to  a  unique  solution  with  a  guarantee  of 
mean  squared-monotonic  convergence  [10-12],  It  is  extremely 
important  to  note  that  while  the  four  sub-matrices  in  [1]  were 
originally  defined  in  terms  of  the  system’s  excitation 
(illuminated  and  shadowed  sides),  the  contraction  mapping 
characteristics  of  the  system  are  independent  of  the  system’s 
excitation.  All  that  is  important  is  that  the  matrix  is  partitioned 
such  that  the  largest  eigenvalue  of  each  sub-matrix  is  less  than 
one.  If  this  case  is  true,  then  by  the  definition  of  a  contraction 
mapping,  the  system  of  linear  equations  may  be  iterated  directly 
to  a  unique  solution  with  a  guarantee  of  mean  squared- 

monotonic  convergence . 

Maue’s  Formulation  LI  1 


Circular  Cylinder  -  Radius  0.5 


Maue’s  Formulation  L12 

Circular  Cylinder  -  Radius  0.5 
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Maue’s  Formulation  L21 


Circular  Cylinder  -  Radius  0.5 


(C) 


Maue’s  Formulation  L22 

Circular  Cylinder  -  Radius  0.5 


Figure  4  -  Eigenvalue  Analysis  the  Hybrid  Iterative  Method 

While  the  guaranty  of  convergence  is  independent  of  excitation 
function  and  the  initial  guess  used  to  start  the  iteration 
‘procedure,  the  rate  of  convergence  is  not.  In  fact  the  rate  of 
convergence  is  heavily  dependent  upon  the  initial  guess  used  to 
start  the  iteration  procedure.  To  illustrate  this  point,  the  current 
distribution  on  a  0.5  X  radius  circular  cylinder  is  solved  for  an 
incident  plane  wave  using  two  different  initial  guesses.  For  this 
particular  example,  the  Total  Field  formulation  was  utilized 
directly  since  all  of  its  eigenvalues  are  less  than  one.  The  first 
initial  guess  used  was  the  traditional  PO  approximation  used  in 
[1-6],  and  the  second  initial  guess  used  was  the  negated 
reflection  of  the  PO  approximation.  The  PO  approximation  is 
considered  to  be  a  very  good  initial  guess,  while  the  negated 
reflection  of  the  PO  approximation  is  considered  to  be  a  very 
poor  initial  guess.  In  both  cases,  the  exact  solution  was 


obtained  directly  using  MoM.  Figures  5a  and  5b  show  the 
converging  currents  for  each  case. 

Convergence  for  PO  Approximation 


Circular  Cylinder  -  Radius  0.5 


(a) 


Convergence  for  Inverted  PO  Approximation 


Circular  Cylinder  -  Radius  0.5 


(b) 

Figure  5  -  Effects  of  Initial  Guess  on  Convergence  Rate 


Figure  5a  shows  these  currents  for  an  initial  guess  that  is  the 
usual  physical  optics  estimate  of  the  current;  twice  the  incident 
field  on  the  illuminated  side  and  zero  on  the  shadow  side.  Note 
the  vast  improvement  in  the  solution  after  one  iteration,  and  the 
nearly  converged  result  after  five  iterations.  Figure  5b  shows 
the  currents  for  an  initial  guess  of  zero  on  the  illuminated  side 
and  minus  twice  the  incident  field  on  the  shadowed  side;  the 
"  opposite"  of  the  physical  optics  estimate.  Note  while  the 
contraction  mapping  nature  of  the  iteration  procedure 
overcomes  the  poor  initial  guess,  the  rate  of  convergence  is 
significantly  slower,  approximately  five  times  slower  than  for 
the  previous  guess. 
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5.  Summary  and  Conclusions  REFERENCES 


Maue’s  formulation  and  the  Total  Field  formulation  of  the  [1] 
MFIE  have  been  analyzed  as  to  their  suitability  for  iterative 
solutions.  It  was  shown  that  both  formulations  could  be 
direcdy  iterated  to  a  unique  solution  with  a  guarantee  of  mean  [2] 
squared-monotonic  convergence  provided  the  eigenvalues  of 
their  corresponding  matrix  expansions  were  less  than  one.  This 
convergence  is  independent  of  the  initial  guess  used  to  initiate 
the  iterative  process.  The  guarantee  of  monotonic  mean  square 
convergence  is  a  direct  result  of  the  contraction  mapping  [3] 
properties  of  Maue’s  formulation  and  Total  Field  formulation. 

The  conditions  under  which  these  formulations  have 
contraction  mapping  properties  were  developed  in  Section  2  for 
the  integral  forms  and  in  Section  3  for  the  discrete  forms. 

Section  3  showed  how  the  eigenvalues  of  the  matrix 
representation  could  be  used  to  determine  the  convergence 
characteristics  of  a  matrix  used  in  the  iterative  solution  process. 

For  the  cases  presented,  the  eigenvalue  plots  showed  a 
preference  for  the  Total  Field  formulation  over  Maue’s 
formulation  since  the  latter  tended  to  have  larger  eigenvalues 
than  the  former.  Eigenvalues  of  value  greater  than  unity 
indicate  that  the  solution  will  not  have  contraction  mapping 
properties.  Eigenvalues  of  value  equal  to  one  correspond  to 
internally  resonant  cases.  Eigenvalues  of  value  equal  to  zero 
can  allow  the  existence  of  spurious  modes.  In  Section  4,  it  was 
noted  that  the  previously  published  Hybrid  Iterative  Method  [1] 
converged  because  the  division  of  the  geometry  into  two  parts 
produced  (for  the  cases  examined)  sub-matrices  whose 
eigenvalues  were  less  than  unity,  thereby  creating  a  series  of 
contraction  mappings.  The  authors  believe  that  the 

convergence  obtained  in  the  other  direct  surface  current- 
iterative  methods  [2-6]  discussed  occurred  because  of  a  similar 
phenomenon. 

Lastly,  the  material  in  this  paper  shows  how  to  use  either 
formulation  of  the  MFIE  successfully  to  obtain  an  iterative 
solution.  Also,  iterative  solutions  to  the  MFIE  have  the 
potential  to  substantially  reduce  the  required  computational 
time  for  large  PEC  bodies  provided  that  the  iterative  process  is 
initiated  with  a  “good”  initial  guess,  such  as  the  PO 
approximation.  Finally,  additional  computational  speed  can  be 
obtained  for  extremely  large  PEC  bodies  by  the  direct 
parallelization  of  the  iterative  methods  discussed  in  Sections  3 
and  4. 
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Modeling  of  Low-Gain  Antennas  on  Aircraft  Using  APATCH 
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ABSTRACT 

Radiation  patterns  for  low-gain  antennas  such  as 
those  used  for  telemetry  and  collision  avoidance  systems 
were  computed  using  the  APATCH  program.  APATCH 
uses  shooting  and  bouncing  rays  (SBR)  to  compute  the 
radiation  pattern  for  the  antenna  installed  on  a  scattering 
geometry.  A  built-in  antenna  model  was  used  to 
represent  a  telemetry  antenna  and  its  accuracy  verified 
by  comparison  with  measurements  and  results  from  a 
method  of  moments  patch  code.  Antenna  patterns  were 
computed  for  various  locations  on  a  Cessna  172  and  an 
F-18  Hornet. 

I.  INTRODUCTION 

When  an  antenna  is  installed  on  an  aircraft,  its 
radiation  pattern  will  change  due  to  the  interaction 
between  the  radiating  element  and  the  aircraft  surface. 
Scattering  mechanisms,  such  as  single  and  multiple 
reflections,  surface  waves,  and  edge  diffraction,  are 
responsible  for  altering  the  overall  radiation  pattern  of 
the  installed  antenna  [1].  The  installation  of  new 
systems  or  physically  reconfiguring  an  aircraft  can  affect 
the  pattern  performance  of  antennas  already  on  the 
aircraft.  This  occurs  frequently  in  flight  testing,  when 
instrumentation  systems  are  temporarily  installed  on  an 
aircraft.  A  performance  analysis  must  be  conducted  to 
assure  that  die  antennas  provide  adequate  coverage  for 
the  anticipated  aircraft  maneuvers.  Lost  data  may  cause 
a  test  to  be  rescheduled,  leading  to  cost  overruns  and 
schedule  delays. 

This  paper  demonstrates  that  a  simple  current 
source  is  an  accurate  model  for  low-gain  wire  antennas 
of  the  type  commonly  used  for  telemetry.  The  ensuing 
analysis  presents  the  performance  of  two  examples  of 
low-gain  antennas.  The  first  is  a  telemetry  antenna, 
which  is  used  with  an  airborne  instrumentation  system, 
and  the  second  is  an  antenna  being  considered  for  a 
traffic  alert/collision  avoidance  system. 

APATCH,  a  computer  program  developed  by 
DEMACO  Inc.  [2],  uses  shooting  and  bouncing  rays 
(SBR)  to  compute  the  radiation  pattern  of  an  antenna 
installed  on  an  aircraft  or  any  other  electrically-large 
electromagnetic  scatterer  [3].  APATCH  uses  the  same 
SBR  model  as  the  radar  cross  section  prediction  code 
XPATCH,  which  has  been  validated  and  used 
extensively  [4,5,6].  Rays  are  shot  from  the  antenna  with 
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a  weighting  determined  by  its  pattern.  APATCH  has  the 
capability  to  model  several  simple  radiating  elements, 
arrays  of  elements,  or  pattern  data  can  be  provided  in 
tabular  format.  The  Advanced  Computer  Aided  Design 
(ACAD)  program  was  used  to  generate  a  triangular  facet 
model  of  the  aircraft  [7]. 

II.  VALIDATION  OF  THE  ANTENNA  MODEL 

To  evaluate  the  low-gain  antenna  modeled  in 
APATCH,  a  test  case  was  run  for  a  telemetry  antenna 
located  on  a  cylinder.  The  antenna’s  construction  is  a 
quarter-wavelength  wire  that  is  excited  from  one  end. 
The  wire  is  inclined  at  a  45°  angle  with  respect  to  the 
horizontal,  and  the  assembly  is  encased  in  an 
aerodynamically  shaped  radome  as  shown  in  Figure  1. 

Measured  data  for  the  antenna  mounted  on  an  8- 
inch  (20.3  cm)  diameter  by  30-inch  (76.2  cm)  length 
cylinder  was  available  from  Haigh-Farr.  This  data  was 
used  to  evaluate  several  approaches  to  modeling  the 
antenna  in  APATCH.  The  faceted  cylinder  model 
referenced  in  the  APATCH  input  file  is  also  shown  in 
Figure  1.  The  telemetry  antenna  is  situated  halfway 
along  the  length  of  the  cylinder  as  indicated  in  the  figure. 

The  telemetry  antenna  was  modeled  as  a  current 
source,  which  is  one  of  the  built-in  antenna  classes 
available  in  the  APATCH  program.  It  is  a  point  source 
that  launches  rays  with  a  dipole  weighting  [2].  Since 
this  antenna  class  makes  no  provision  for  the  control  of 
the  length  and  radius  of  the  radiating  element,  the  height 
of  the  current  source  ( h )  was  varied  until  a  good  match 
was  achieved  between  the  computed  results  and  the 
measured  data.  This  provides  a  limited  means  by  which 
the  effects  of  element  length,  radius,  and  radome 
properties  can  be  included. 

Figure  2  compares  measured  data  with  the 
APATCH  results  for  two  heights.  The  agreement 
between  measured  and  computed  data  is  generally  better 
in  both  the  upper  and  lower  hemispheres  for  a  height  of 
h-\l  mm  as  compared  to  a  height  of  h=2  mm. 

The  x-y  plane  pattern  computed  by  APATCH 
closely  resembles  the  measured  data  as  shown  in  Figure 
2b.  The  two  curves  have  been  shifted  from  the 
normalized  value  for  clarity.  The  symmetric  behavior  of 
the  scattered  field  about  the  scatterer  is  exhibited  in  these 
plots. 
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Figure  1:  Telemetry  antenna  and  ACAD  model  of  cylinder  used  in  the  preliminary  study. 


The  field  due  to  edge  diffraction  can  be  computed 
by  APATCH  if  an  edge  file  is  provided.  The  edge  file  is 
a  sequential  list  containing  edge  length,  orientation,  and 
wedge  angle  for  all  edges  in  the  ACAD  model.  This 
APATCH  utility  was  implemented,  but  it  did  not 
significantly  improve  the  agreement  between  the 
computed  and  measured  results.  This  may  have  been 
due  to  the  relatively  coarse  mesh  of  the  cylinder. 

The  method  of  moments  code  PATCH  [8]  was  also 
used  to  verify  the  antenna  and  cylinder  model.  PATCH 
uses  the  triangular  surface  subdomains  with  overlapping 
rooftop  basis  functions  [9].  Unlike  APATCH,  where  the 
antenna  pattern  is  decoupled  from  the  scatterer,  PATCH 
includes  interactions  between  the  cylinder  and  antenna. 
Figure  3  shows  a  detail  of  the  region  where  the  antenna 
is  attached  to  the  cylinder  in  the  ACAD  model.  One  of 
the  edges  of  the  facets  comprising  the  antenna  is  excited 


with  a  voltage  source,  and  PATCH  computes  the 
resulting  surface  currents  and  scattered  field. 

Figure  3  also  shows  the  elevation  plane  pattern. 
PATCH  computed  the  gain  to  be  4.1  dB  (compared  with 
4.9  dB  from  APATCH).  Because  PATCH  is  a  method- 
of-moments  code,  from  a  practical  perspective,  it  cannot 
be  used  to  compute  the  radiation  patterns  of  antennas  on 
electrically  large  aircraft.  Computation  times  for 
problems  involving  large  scattering  geometries  can  be 
quite  high.  Furthermore,  the  demand  for  computer 
memory  rises  very  quickly  as  the  electrical  size  of  the 
scattering  geometry  increases.  The  close  agreement 
between  PATCH,  APATCH,  and  the  measured  data 
indicates  that  the  APATCH  current  source  antenna 
model  is  a  good  representation  of  the  actual  stub. 


(a) 


(b) 


Figure  2:  APATCH  results  compared  with  measured  data  for  the  cylinder  and  telemetry  antenna.  Ee  polarization  shown, 
(a)  y-z  plane  (b)  x-y  plane 
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Figure  3:  PATCH  model  and  results  for  cylinder  with  telemetry  antenna,  (a)  detail  of  ACAD  model  used  in  PATCH 
computation,  (b)  elevation  plane  pattern 


Figure  4:  ACAD  model  of  F-18  Hornet  (2620  facets).  Locations  where  the  telemetry  antenna  was  considered  are 
designated  as  A  and  B. 


IIL  F-18  HORNET  WITH  TELEMETRY 
ANTENNA 

The  antenna  model  discussed  in  Section  II  was 
simulated  installed  on  the  F-18  aircraft  shown  in  Figure 
4.  With  a  highly  dynamic  aircraft  such  as  the  F-18,  it  is 
likely  that  two  antennas  would  be  used  in  a  telemetry 
link:  one  on  the  bottom  and  one  on  top.  This  approach 
ensures  that  there  is  adequate  reception  at  the  test 


mission  control  facility  as  the  aircraft  executes  the 
maneuvers  required  for  a  particular  test.  The  two 
antennas  are  connected  to  the  transmitter  by  a 
microwave  coupler.  Usually  a  3  dB  coupler  is  used 
which  equally  divides  the  signal  power  into  the  two 
antennas.  However,  there  are  situations  where  the 
majority  of  the  time  the  lower  antenna  is  visible  to  the 
control  facility.  In  these  cases  it  may  be  more  efficient 
to  use  a  10  dB  coupler  so  that  10%  of  the  input  power  is 
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directed  to  the  upper  antenna  and  90%  to  the  lower 
antenna-  The  upper  antenna  is  supplied  so  that  the 
telemetry  link  is  maintained  during  the  brief  time  that  a 
maneuver  is  executed,  albeit  the  power  from  this  antenna 
is  significantly  reduced. 

APATCH  results  are  shown  in  Figure  5  for  the  10 
dB  coupler.  The  Ee  and  E $  polarizations  are  shown  for 
the  roll,  pitch,  and  yaw  planes.  The  reference  gain 
values  (i.e.,  the  0  dB  levels)  are  listed  for  each  plot. 
(Note  that  insignificant  cross-polarized  values  were  not 
plotted.  For  example,  there  is  no  significant  E $  in  the 
pitch  plane.)  It  can  be  seen  that  there  is  good  coverage 
in  the  lower  region.  The  region  above  the  aircraft 
exhibits  lower  radiation  levels  than  the  region  below  the 
aircraft;  this  is  expected  when  the  10  dB  coupler  is  used. 
Finally,  the  yaw  plane  pattern  is  shown  in  Figure  5c.  It 
displays  the  interference  between  the  upper  and  lower 
antennas  characterized  by  its  oscillatory  shape. 

When  a  3  dB  coupler  is  used  with  the 
instrumentation  system,  the  effective  radiation  patterns 
from  the  upper  and  lower  telemetry  antennas  are  shown 
in  Figure  6.  The  roll  and  pitch  plane  patterns  indicate 
that  there  is  adequate  radiation  into  the  regions  above 
and  below  the  aircraft.  This  is  the  type  of  coverage 
desired  when  the  aircraft  is  performing  test  missions 
requiring  frequent  maneuvers  and  unusual  attitudes  like 
those  encountered  in  tests  of  performance  and  flying 
qualities. 

The  yaw  plane  pattern  shown  in  Figure  6c  displays 
large  variations  in  the  E&  and  E$  polarizations  in  the 
horizontal  plane.  The  peaks  and  nulls  of  the  two 
polarizations  occur  in  an  alternating  fashion  due  to  the 
interferometer  effect  of  two  widely  separated  elements. 
This  effect  is  more  pronounced  for  the  3  dB  coupler  than 
for  the  10  dB  coupler  because  the  equal  power  weighting 
provides  nearly  complete  cancellation. 

IV.  CESSNA  172  WITH  TCAS  ANTENNA 

An  antenna  suitable  for  use  with  the  Traffic  Alert 
and  Collision  Avoidance  System  (TCAS)  was  simulated 
on  the  Cessna  172  (C-172)  aircraft.  The  US  Federal 
Aviation  Administration  (FAA)  mandates  the  use  of 
TCAS  on  all  passenger  and  cargo  aircraft  with  10  seats 
or  more.  TCAS  aids  pilots  in  detecting  the  presence  of 
nearby  aircraft.  The  system  operates  by  interrogating  the 
transponder  of  an  approaching  aircraft.  Range  is 
determined  from  the  time  required  to  receive  the  reply, 
and  altitude  is  extracted  from  the  information  encoded 
within  the  reply  [10].  Mode  C  is  the  term  used  to 
describe  a  transponder  that  reports  its  altitude  when 
interrogated.  In  a  complete  TCAS  system,  a  pair  of 
directional  antennas  is  used  to  establish  the  bearing  of 


the  approaching  aircraft.  This  version  is  not  explored 
here.  Instead,  only  an  omni-directional  antenna  is 
considered,  which  only  yields  the  range  and  altitude  of 
an  intruding  aircraft.  It  is  proposed  that  bearing  be 
determined  from  a  position  reporting  feature  of  a 
transponder  similar  to  Mode  C  operation.  TCAS 
operates  on  two  frequencies.  It  interrogates  (transmits) 
on  1030  MHz  and  receives  transponder  replies  on  1090 
MHz. 

The  aircraft  model  used  in  the  following  analysis  is 
shown  in  Figure  7.  Indicated  in  the  figure  is  the  location 
considered  for  antenna  placement.  A  vertically  polarized 
quarter-wavelength  monopole  antenna  was  used,  which 
in  free  space  produces  an  omni-directional  radiation 
pattern.  Figure  8  shows  the  radiation  pattern  for  the 
quarter-wave  element  installed  on  the  aircraft.  The  roll 
pattern,  Figure  8a,  shows  good  coverage  in  the  region 
above  the  aircraft,  except  for  directly  above,  where  there 
exists  a  sharp  null  in  the  pattern  which  is  typical  of 
vertically  polarized  antennas.  There  is  less  coverage  in 
the  region  from  the  horizon  to  approximately  30°  below 
the  horizon.  This  is  the  effect  of  the  wings  on  the  pattern 
since  the  antenna  is  masked  by  the  large  wing  surfaces. 

Referring  to  the  pitch  plane  pattern.  Figure  8b,  the 
gain  in  the  aft  region  is  higher  than  that  in  the  front.  The 
yaw  plane  pattern  shows  uniform  radiation  in  the 
horizontal  plane  surrounding  the  aircraft.  Antenna  gain 
in  the  forward  direction  is  comparable  to  that  found  in 
the  aft  region.  With  the  antenna  mounted  at  this 
location,  the  interaction  with  the  wings  gives  an 
increased  gain  that  improves  the  detection  range  in  the 
region  above  the  wings  where  the  pilot  experiences 
obstructed  visibility. 

By  incorporating  an  array  antenna,  the  gain  in  the 
aft  region  can  be  reduced  in  exchange  for  greater 
coverage  in  the  forward  region;  this  would  improve 
TCAS  performance  and  aircraft  safety  even  more.  A 
simple  array  can  be  manufactured  consisting  of  two 
vertical  elements,  each  a  quarter- wavelength  long  and 
separated  by  one-quarter  wavelength.  The  aft  element  is 
fed  90°  out  of  phase  with  respect  to  the  forward  element 
to  obtain  a  cardioid  shape.  Figure  9  shows  the 
performance  of  the  array  antenna  installed  on  this 
aircraft. 

The  roll  plane  pattern  (Figure  9a)  shows  that  the 
coverage  in  the  space  above  the  aircraft  is  similar  to  the 
single  antenna  at  the  same  location,  but  the  gain  has 
increased  slightly  from  8.14  dB  to  8.56  dB  for  1030 
MHz. 

When  the  plots  of  the  pitch  and  yaw  planes  are 
examined,  the  benefit  of  the  array  becomes  apparent 
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Figure  7:  ACAD  model  of  the  C-172  aircraft  (2790 
facets). 

(see  Figures  9b  and  9c).  An  increase  in  gain  of 
approximately  3  dB  in  the  forward  direction  is  evident. 
Gain  in  this  direction  exceeds  the  gain  in  the  aft  region. 
Additionally,  when  compared  with  the  single  element 
antenna,  a  net  gain  of  0.8  dB  is  achieved  in  the 
horizontal  plane. 

IV.  SUMMARY  AND  CONCLUSIONS 

APATCH  is  a  code  used  to  compute  the  radiation 
patterns  of  antennas  installed  on  complex  objects.  It  is 
based  on  the  SBR  technique,  which  is  applicable  to 
electrically  large  structures.  The  built-in  current  element 
was  used  to  model  an  antenna  on  a  cylindrical  test 
fixture  and  the  C-172  and  F-18  aircraft.  Comparison  of 
the  APATCH  results  with  measured  data  and  the  results 
from  a  method  of  moments  patch  code  shows  that  the 
built-in  current  source  is  an  accurate  model  for  low-gain 
wire  antennas  of  the  type  commonly  used  for  telemetry. 
The  run  times  for  the  F-18  and  C-172  were 
approximately  17  and  5  minutes,  respectively,  on  a  SGI 
Indigo2.  CPU  times  can  vary  widely  depending  on  the 
number  of  facets,  ray  density,  and  number  of  bounces. 

On  a  maneuvering  F-18,  good  telemetry  link 
performance  can  be  expected  from  a  pair  of  antennas 
mounted  on  the  top  and  bottom  of  the  aircraft.  When  a  3 
dB  coupler  is  used  to  split  the  telemetry  signal,  an 
interferometer  pattern  results  in  the  yaw  plane.  This 
behavior  can  lead  to  telemetry  drop-outs  when  a  linearly 
polarized  antenna  is  used  in  the  receiving  ground  station. 

This  paper  also  presented  a  performance  analysis  of 
an  antenna  considered  for  a  type  of  TCAS  system  on  the 
C-172  aircraft.  Acceptable  performance  can  be  achieved 
when  either  a  single  quarter-wave  element  or  a  two- 
element  array  is  used. 
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Figure  8:  C-172  aircraft  with  TCAS  antenna 
(a)  roll,  (b)  pitch,  and  (c)  yaw 


.7\> 


* 


—  fo-1 030  MHz  04.07  dB  270 

-.f  -1090  MHz  G-4.02  dB 


Figure  9:  C-172  aircraft  with  array  antenna 
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Abstract — This  paper  describes  a  new  contribution 
to  the  analysis  of  arbitrary  shielded  circuits  and  an¬ 
tennas  of  complex  shapes,  in  the  frame  of  the  integral 
equation  (EE)  and  Method  of  Moments  formulation 
(MoM).  The  technique  is  based  on  the  spatial  image 
approach  and  a  new  specially  truncated  set  of  images 
is  developed  to  enhance  the  convergence  behavior  of 
the  series  involved.  Results  show  that,  with  the  new 
specially  truncated  series  of  images,  convergence  is 
achieved  very  fast.  In  this  paper  simulated  results 
obtained  with  the  new  approach  are  compared  with 
measurements. 

Index  terms — Microwave  circuits  and  antennas, 
shielding,  Green  function,  spatial  images,  integral 
equation,  method  of  moments. 

I.  Introduction 

The  analysis  of  shielded  microwave  circuits  and  anten¬ 
nas  is  a  subject  that  has  always  attracted  much  attention 
and  several  numerical  models  have  been  developed  in  the 
past.  Among  them,  finite  element  techniques  have  been 
successfully  used  [1]  but  they  usually  lead  to  computer 
codes  which  are  computationally  heavy.  Perhaps  the  most 
popular  and  efficient  technique  is  the  integral  equation 
(IE)  formulation  combined  with  the  Method  of  Moments 
(MoM)  algorithm  [2].  The  main  practical  difficulty  in  this 
approach,  however,  is  in  the  numerical  evaluation  of  the 
associated  Green’s  functions  usually  formulated  as  very 
slowly  convergent  modal  series.  For  the  efficient  summa¬ 
tion  of  these  series  the  Fast  Fourier  Transform  (FFT)  has 
been  successfully  used  in  the  past  [2], [3],  but  it  restricts 
the  subsequent  discretization  of  the  circuits  to  uniform 
meshes.  Other  acceleration  techniques,  without  the  use 
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of  the  FFT,  have  been  recently  reported,  namely  the  use 
of  the  residue  theorem  [4]  and  the  use  of  the  Sommerfeld 
identity  [5].  An  alternative  technique  to  the  modal  series 
approach  is  to  formulate  the  Green’s  functions  in  the  spa¬ 
tial  domain  using  the  image  theory.  In  this  context  we 
can  mention  the  work  in  [6]  who  used  the  Ewalt  transfor¬ 
mation  to  sum  the  spatial  images  series,  but  no  dielectric 
layers  are  considered,  and  the  work  in  [7]  who  formulated 
the  spatial  images  series  only  to  the  quasi-static  part  of 
the  kernel. 

In  the  frame  of  IE- MoM,  we  provide  in  this  paper  a 
new  contribution  to  the  analysis  of  arbitrary  shielded  cir¬ 
cuits  and  antennas  of  complex  shapes.  The  technique  is 
based  on  the  spatial  images  approach  and  it  formulates 
the  complete  Green’s  functions  associated  to  a  general 
multilayered  medium  as  standard  Sommerfeld  integrals. 
Once  Sommerfeld  type  Green’s  functions  are  computed, 
the  effect  of  lateral  walls  is  locally 'Included  by  adding 
spatial  images.  Care  is  exercised  to  ascertain  the  conver¬ 
gence  behavior  of  the  resulting  spatial  images  series.  For 
those  cases  exhibiting  slow  convergence  behavior,  a  spe¬ 
cially  truncated  series  of  images  is  developed  based  .on  the 
imposition  of  the  boundary  conditions  for  the  fields  and 
potentials  at  specific  points  on  the  metallic  walls.  Results 
reveal  that  convergence  with  the  specially  truncated  set 
of  images  is  very  fast  and  accurate  results  are  obtained 
with  only  10  spatial  images.  An  interesting  feature  of 
the  technique  developed  is  that  triangular-cell  based  MoM 
formulations  can  be  easily  used,  thus  allowing  the  analy¬ 
sis  of  arbitrary  complex  structures.  In  addition,  with  the 
proposed  approach  all  the  know-how  in  the  field  of  Som¬ 
merfeld  integrals  evaluation  can  be  directly  put  to  work 
for  the  analysis  of  shielded  structures. 

II.  Image  theory 

The  key  step  in  the  analysis  of  planar  multilayered 
printed  circuits  following  the  space  IE  approach,  is  the 
derivation  of  the  spatial  domain  Green’s  functions  formu¬ 
lated,  in  the  case  of  a  medium  of  infinite  lateral  trans¬ 
verse  dimensions,  as  Sommerfeld  integrals  [8].  With  ref- 


1054-4887  ©  1999  ACES 


92 


ACES  JOURNAL,  VOL  14,  NO.  3,  NOVEMBER  1999 


erence  to  the  shielded  structure  in  Fig.  1,  the  correspond¬ 
ing  Green’s  functions  accounting  for  the  presence  of  the 
lateral  walls  can  be  written  using  standard  spatial  image 
theory.  The  Green’s  functions  are  therefore  again  formu- 
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Fig.  1.  Basic  shielded  multilayered  structure  analyzed  in  this  paper 
and  the  associated  spatial  images  for  a  unit  point  charge. 

lated  as  Sommerfeld  integrals,  but  they  are  computed  for 
an  extended  range  of  source-observer  distances  to  account 
for  the  interaction  of  all  required  spatial  images.  For  this 
purpose  asymptotic  techniques  specially  efficient  for  large 
source-observer  distances  are  best  utilized  as  described  in 
[9].  The  boxed  mixed  potential  Green’s  functions  GBox 
are  then  written  as  [8] 

GBOx(*Ej  y\x  j  y )  — 

oo  oo 

y;  y  \G(x,y\+x'  +  2ma,+y'  +  2nb) 

m— — oo  n— — oo 

+sx  G(x,  y\  -x'  +  2ma,  +  yf  4-  2nb) 

+sy  G(x ,  y\  +x'  4-  2  ma,  -y'  +  2  nb) 

sy  G(x ,  y\  -x'  4-  2 ma,  -yf  4-  2n6)j , 

where  G  is  the  corresponding  Sommerfeld  type  Green’s 
function  and  sy  are  sign  functions  taking  the  values 
shown  in  Table  I  for  the  different  types  of  mixed  poten¬ 
tial  Green’s  functions  components.  As  shown  in  (1)  the 


total  Green’s  functions  are  now  expressed  in  terms  of  in¬ 
finite  series  and  the  convergence  behavior  becomes  there¬ 
fore  an  issue.  For  instance,  we  have  found  convergence 
problems  in  the  analysis  of  the  microstrip  line  shown  in 
Fig.  2,  printed  on  a  thick  dielectric  substrate  (er  =  3.0, 
h/X  =  0.2). 

To  illustrate  this  point,  we  present  in  Fig.  3(a)  and 
3(b)  the  computed  S-parameters  when,  respectively,  22 
and  102  spatial  images  are  included  in  the  analysis.  As 
it  can  be  observed,  the  ripple  in  the  computed  response 
indicates  that  convergence  rapidly  deteriorates  with  fre¬ 
quency.  Moreover,  the  use  of  more  images  does  not  solve 
the  problem  and,  as  shown  in  Fig.  3(b),  only  the  ripple 
is  more  compressed.  This  slowly  convergence  behavior 
can  be  explained  as  being  due  to  the  excitation  of  surface 
waves  in  the  structure.  To  easily  understand  this  fact  we 
have  computed  the  associated  Green’s  functions  using  the 
imaginary  axis  decomposition  reported  in  [9].  Following 
this  technique,  the  total  Green’s  functions  are  expressed 
as  the  sum  of  three  main  contributions,  namely  the  quasi¬ 
static  part,  the  spatial  wave  part  and  the  contribution  due 
to  the  surface  waves  excitation.  Fig.  4  shows  the  com¬ 
puted  electric  scalar  potential  Green’s  function  for  two 
different  frequencies.  The  first  one  (2  GHz)  corresponds 
to  a  frequency  where  convergence  is  still  good  while  the 
second  one  (10  GHz)  corresponds  to  a  frequency  well  in¬ 
side  the  ripple  region.  We  can  clearly  observe  that  at  2 
GHz  the  surface  wave  is  very  weakly  excited  and  it  starts 
to  dominate  the  global  Green’s  function  behavior  for  spa¬ 
tial  distances  (p)  of  order  of  ( kop  =  100).  Moreover, 
for  small  source-observer  distances  the  Green’s  function 
exhibit  a  strong  decaying  behavior  of  (1/p2)  type.  On 
the  contrary,  at  10  GHz  the  surface  wave  is  very  strongly 
excited  since  it  dominates  the  global  Green’s  function  be¬ 
havior  for  spatial  distances  of  about  only  (ho  p  =  1).  In 
addition,  for  small  source-observer  distances  the  decaying 
behavior  is  now  inverse  of  the  distance  (1/p),  instead  of 
(1/p2)  as  before.  Consequently,  we  can  conclude  that  the 
convergence  properties  of  the  spatial  images  series  are  re¬ 
lated  to  the  amount  of  electromagnetic  energy  coupled  to 
the  surface  wave  modes.  To  try  to  overcome  this  problem, 
a  specially  truncated  set  of  images  has  been  developed  and 
it  is  described  in  the  next  section. 

TABLE  I 

Value  of  the  sign  functions  for  all  different  mixed  potential  Green’s  functions 
components.  G-qA  denotes  boxed  magnetic  vector  potential  Green’s  function 
while  G-qv  is  the  boxed  electric  scalar  potential. 
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III.  Specially  Truncated  Set  of  Images 

The  problem  of  the  slow  convergence  behavior  of  the 
series  developed  in  the  previous  section  can  be  overcome 
by  noticing  that  there  are  two  main  features  in  any  boxed 
Green’s  function  that  must  be  preserved  in  order  to  obtain 
accurate  results,  namely: 

1.  The  singular  behavior  when  p  — >  0. 

2.  The  boundary  conditions  at  all  lateral  cavity  walls. 

The  singular  behavior  is  naturally  preserved  in  the  devel¬ 
oped  series  of  images,  because  they  have  been  constructed 
using  standard  Green’s  functions.  What  remains  then, 
is  the  accurate  imposition  of  the  boundary  conditions  at 
the  metallic  walls.  This  can  be  done  by  simply  adding 
the  first  few  images  of  the  series  as  before,  and  then  com¬ 
puting  the  needed  strength  of  the  last  remaining  images 
so  that  the  boundary  conditions  for  the  fields  are  strictly 
satisfied  at  the  projections  of  the  observer  points  at  the 
metallic  walls.  To  simplify  the  analytical  exposition  of  the 
problem  we  now  show  the  procedure  for  the  case  of  the 
electric  scalar  potential  inside  a  parallel  plate  waveguide, 
being  formally  analogous  for  the  rectangular  cavity  case. 
Prom  (1)  and  Table  I  we  write  the  electric  scalar  potential 
Green’s  function  in  the  parallel  plate  waveguide  structure 
shown  in  Fig.  5,  as 

GBV(x,y\x',y')  = 

OO 

^2  [Gv(x,y\+x'+2ma,y')  (2) 

m=- oo 

-Gv(x,y\  -x'+2ma,2//)]J 

where  Gv  is  the  Sommerfeld  type  electric  scalar  potential 
Green’s  function  in  the  layered  structure  considered.  In 
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(a)  Results  with  22  spatial  images. 


(b)  Results  with  102  spatial  images. 

Fig.  3.  Scattering  parameters  of  the  structure  in  Fig.  2  when  22  and 
102  spatial  images  are  included  in  the  analysis 


Fig.  2.  Microstrip  line  printed  on  a  thick  substrate  and  backed  by 
a  parallel  plate  waveguide. 
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the  above  equation,  it  is  convenient  to  define  an  interme¬ 
diate  Green’s  function  representing  the  potential  created 
by  a  basic  image  couple  (basic  image  set)  as  (Fig.  5) 
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(a)  Frequency  2  GHz. 


(b)  Frequency  10  GHz. 

Fig.  4.  Imaginary  axis  decomposition  for  the  electric  scalar  potential 
associated  to  the  problem  in  Fig.  2.  Source-observer  distance  is  p. 


GV*(x,y\x',y')  =Gv(x,y\  +x'+2ma,y')  - 
Gv(x,y\  -x'+2ma,y'), 

and  obviously,  this  intermediate  Green’s  function  matches 
the  boundary  conditions  at  the  x  —  0  wall.  To  serve  our 
purposes,  it  is  convenient  to  truncate  the  infinite  spatial 
series  with  (M)  and  rewrite  equation  (2)  using  the  defini¬ 
tion  in  (3)  as 

+M-1 

GBV(x,y\x',y')&  ^  G^{x,y\x',y') 

771—  —  M  -)-l 

+  G%M  (x,  y\x’,  y’)  +  GV!m  (x,  y\x',  y /).  (4) 

It  is  interesting  to  observe  that  the  spatial  images  in  (4) 
is  balanced  with  respect  the  first  metallic  wall  at  (x  =  0), 
that  is  the  boundary  conditions  are  already  rigorously  im¬ 
posed  at  this  wall.  On  the  contrary,  the  series  is  unbal¬ 
anced  with  respect  the  second  metallic  wall  at  (x  =  a). 
The  key  idea  in  this  technique  is  to  introduce  an  unknown 
constant  ( q )  so  that  the  strength  of  the  two  last  basic  im¬ 
age  sets  in  the  series  are  adjusted  to  enforce  the  right 
boundary  condition  for  the  potential  also  at  the  second 
metallic  wall.  We  then  introduce  the  unknown  constant 
(i q )  and  rewrite  (4)  as 

+M-1 

GBV(x,y\x',y')&  G$*(x,y\x',y') 

771=—  M+l 

+  q  (x,  y\x',  y')  +  G^fM  (x,  y\x’,  y')]  •  (5) 


Fig.  5.  Spatial  images  arrangement  for  a  parallel  plate  waveguide 
structure. 
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In  (5)  the  potential  created  by  the  first  (2 M  —  1)  basic 
image  sets  can  now  be  defined  as 


+M- 1 

GBV{x,y\x',y')  =  ^  Gt$*(x,y\x',  y'),  (6) 

and  the  potential  created  by  the  two  last  basic  image  sets 
as 

GlBV(x,y\x',y')  =  q  [Gv*M(x,y\x’,y')  +  G$*M(x,y\x' ,y')] . 

(7) 

Next,  the  potential  created  by  the  first  (2 M  -  1)  basic 
image  sets  at  the  projection  of  the  observer  point  at  the 
second  metallic  wall  is  computed  from  (6)  as 

+M- 1 

w  =  GBV(a,y\x’ ,y')  =  G^*(a,y\x',y').  (8) 

m——M 4*1 

Since  the  total  potential  must  be  zero  also  at  the  second 
metallic  wall,  the  potential  created  by  the  two  last  basic 
image  sets  must  compensate  (8),  and  thus  the  following 
relation  can  be  written 

-w  =  q  [g®^ (a, y\x\ y')  +  G^M (a, y\x’, y')] ,  (9) 


(a)  Frequency  10  GHz. 


from  where  the  unknown  constant  required  to  annihilate 
the  total  potential  also  at  the  second  metallic  wall  is  easily 
found  as 


—w 

7i  +  72  ’ 


(10) 


with  the  redefinitions  of  the  following  coefficients 


7i  =  G^sM(a,y|x',y,); 


72  =  Gf*M(a,y\x',y').  (11) 


It  is  important  to  point  out  that  a  similar  procedure 
as  the  one  described  here  for  the  electric  scalar  potential 
is  also  applied  for  those  components  of  the  magnetic  vec¬ 
tor  potential  satisfying  Newman  boundary  conditions  at 
the  metallic  walls,  and  the  only  difference  is  that  numer¬ 
ical  differentiation  is  used  in  order  to  annihilate  the  total 
derivative. 

To  show  the  effectiveness  of  the  approach  developed,  we 
have  computed  the  electric  scalar  potential  Green’s  func¬ 
tion  for  the  structure  in  Fig.  2  when  the  direct  sum  in 
(4)  is  used  and  when  the  new  specially  truncated  set  of 
images  is  used.  In  Fig.  6  we  show  the  results  obtained 
as  a  function  of  the  observer  point  position  for  two  dif¬ 
ferent  frequencies  (10  GHz  and  30  GHz),  fixed  position  of 
the  source  (x;/a  =  0.107)  and  when  a  total  number  of  3 
and  60  pairs  of  basic  image  sets  are  included  in  the  direct 
sum,  and  only  3  pairs  of  basic  image  sets  in  the  specially 


(b)  Frequency  30  GHz. 


Fig.  6.  Shielded  electric  scalar  potential  Green’s  function  associated 
to  the  problem  in  Fig.  2  when  direct  sum  is  used  and  when  the  new 
specially  truncated  set  of  images  is  used. 
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truncated  series.  As  we  can  observe,  close  to  singularity 
both  approaches  give  similar  results,  but  boundary  condi¬ 
tions  at  the  second  lateral  wall  are  only  satisfied  when  the 
specially  truncated  set  of  images  is  used.  Also  we  observe 
that  the  oscillations  due  to  the  standing  wave  created  by 
the  presence  of  the  lateral  walls  are  slightly  readjusted 
when  the  specially  truncated  set  of  images  is  used.  This 
readjustment  appears  to  be  necessary  to  allow  the  poten¬ 
tial  to  be  zero  at  the  lateral  metallic  walls.  Finally,  it 
can  be  observed  that  when  60  pairs  of  basic  image  sets 
(242  images)  are  used  in  the  direct  sum,  the  boxed  elec¬ 
tric  scalar  potential  converges  to  the  results  obtained  with 
just  3  specially  truncated  pairs  of  basic  image  sets  (14  im¬ 
ages),  thus  indicating  the  huge  saving  in  computational 
effort  obtained  with  the  new  approach. 

When  the  newly  developed  Green’s  functions  are  used 
in  the  analysis  of  the  structure  in  Fig.  2  we  observe  fast 
convergence  properties  in  the  results.  In  Fig.  7  we  present 
the  computed  S-parameters  showing  very  smooth  behav¬ 
ior  and  that  the  ripple  has  been  effectively  canceled  out, 
even  with  only  3  specially  truncated  pairs  of  basic  image 
sets  (14  images).  By  contrast,  convergence  is  not  achieved 
even  after  25  pairs  of  basic  image  sets  (102  images)  with¬ 
out  truncation  as  shown  in  Fig.  3. 

To  further  validate  the  theory  presented,  a  breadboard 
of  the  structure  in  Fig.  2  has  been  manufactured  and 
tested.  In  Fig.  8  we  present  measured  versus  simulated 
results  showing  good  agreement  for  both  Sn  and  S21 
scattering  parameters  and  with  only  3  specially  truncated 
pairs  of  basic  image  sets  in  the  calculations. 

Finally,  it  is  important  to  point  out  that  with  the  pro¬ 
posed  approach  all  additional  computations  are  performed 
analytically  and  that  the  reduction  in  the  number  of  im¬ 
ages  required  to  obtain  good  convergence  is  considerable, 
so  that  the  implementation  of  efficient  software  codes  is 
possible. 


IV.  Results  and  Discussions 

A  software  code  based  on  the  proposed  approach  has 
been  written  for  the  analysis  of  complex  cavity  backed 
antennas.  The  asymptotic  technique  described  in  [9]  has 
been  used  for  the  evaluation  of  the  basic  Sommerfeld  in¬ 
tegrals  for  large  values  of  source-observer  distances  and  a 
MoM  technique  based  on  triangular  cells  has  been  imple¬ 
mented. 

To  show  the  usefulness  of  the  technique,  we  present 
the  analysis  of  a  circular-polarized  cavity-backed  antenna 
containing  two  stacked  patches  of  complex  shapes  (note 
in  particular  the  narrow  slits  existing  in  both  patches).  In 


(a)  Results  with  2  truncated  pairs  of  basic  image  sets. 
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(b)  Results  with  3  truncated  pairs  of  basic  image  sets. 


Fig.  7.  Scattering  parameters  of  the  structure  in  Fig.  2  when  2  and 
3  specially  truncated  pairs  of  basic  image  sets  are  included  in  the 
analysis. 
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(a)  Sn  parameter. 


(b)  521  parameter. 

Fig.  8.  Measured  versus  simulated  results  for  the  scattering  parame¬ 
ters  of  the  structure  in  Fig.  2.  Only  3  truncated  pairs  of  basic  image 
sets  (14  images)  are  used  in  the  theoretical  prediction. 


Fig.  9  we  show  the  basic  geometry  of  the  antenna  and  the 
triangular  meshes  used  in  conjunction  with  the  developed 
algorithms.  For  the  analysis  of  this  structure,  adequate 
magnetic  currents  are  defined  at  the  cavity’s  top  aper¬ 
ture  and  a  set  of  coupled  integral  equations  is  defined  and 
solved  in  the  fashion  described  in  [10].  In  these  equa¬ 
tions,  all  electromagnetic  interactions  are  computed  with 
the  specially  truncated  set  of  images  derived  in  this  pa¬ 
per.  Fig.  10  presents  measured  versus  simulated  results 
for  the  S-parameters  of  the  antenna  when  the  upper  patch 
is  removed,  and  Fig.  11  for  the  complete  antenna  struc¬ 
ture,  with  the  upper  patch  on  top,  showing  in  all  cases 
good  agreement.  An  interesting  feature  of  the  technique 
derived  is  the  rate  of  convergence  of  the  numerical  re¬ 
sults  when  the  number  of  specially  truncated  images  is 
increased.  For  this  last  example  no  significance  numerical 
difference  is  observed  with  10  and  14  images,  thus  indi¬ 
cating  good  and  fast  convergence  of  the  developed  series. 

V.  Conclusions 

In  the  frame  of  IE-MoM,  a  new  contribution  to  the  anal¬ 
ysis  of  arbitrary  shielded  microwave  circuits  and  cavity 
backed  antennas  has  been  presented.  Following  the  pro¬ 
posed  approach,  Green’s  functions  are  formulated  with 
well  known  Sommerfeld  integrals  and  the  know-how  al¬ 
ready  available  in  this  field  is  reused.  The  presence  of 
lateral  metallic  walls  is  taken  into  account  by  using  stan¬ 
dard  image  theory  and  the  convergence  conditions  of  the 
resulting  image  series  have  been  investigated.  When  slow 
convergence  occurs,  a  specially  truncated  set  of  images 
has  been  derived,  based  on  the  rigorous  imposition  of  the 
boundary  conditions  for  the  fields  and  potentials  at  the 
metallic  walls.  Results  reveal  that  convergence  using  the 
specially  truncated  set  of  images  is  very  fast  and  that  with 
only  10-14  images  accurate  results  are  obtained.  More¬ 
over,  the  use  of  the  technique  described  simply  adds  few 
analytical  operations  to  the  basic  Sommerfeld  formalism 
so  that  efficient  computer  codes  can  still  be  developed.  In 
this  paper  measured  and  simulated  results  are  compared 
for  a  stacked  patch  antenna  with  both  patches  exhibiting 
a  complicated  shape. 
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Fig.  11.  Measured  versus  simulated  results  for  the  antenna  in  Fig.  9 
when  the  upper  patch  is  placed  on  top.  Real  and  imaginary  part  of 
the  S- parameters  are  given. 
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Abstract —  A  common  bottle-neck,  limiting  the 
performance  of  many  electromagnetic  numerical 
methods,  is  the  solution  of  sparse  linear  systems. 
Until  now,  this  task  has  been  typically  solved  by 
using  iterative  sparse  solvers,  whose  require  heavy 
computational  efforts,  especially  when  the  problem 
is  not  well- conditioned. 

An  alternative  strategy  is  based  on  the  use 
of  banded  solvers,  which  numerical  complexity  is 
quadratical  with  respect  to  the  matrix  bandwidth. 
Of  course,  these  methods  are  efficient  provided  that 
the  matrix  bandwidth  is  sufficiently  small.  In  this 
paper,  a  method  (called  WBRA)  for  the  bandwidth 
reduction  of  a  sparse  matrix  is  presented:  it  is  here 
specifically  customized  to  typical  electromagnetic 
matrices.  The  approach  is  superior  to  all  the  pre¬ 
vious  algorithms,  also  with  respect  to  commercial 
well-known  packages,  and  is  suitable  also  for  non- 
symmetric  problems. 

As  demonstrated  by  results,  the  use  of  WBRA, 
in  conjunction  with  common  banded  solvers,  sub¬ 
stantially  improves  (up  to  one  order  of  magnitude) 
the  solution  times  in  several  electromagnetic  ap¬ 
proaches,  such  as  Mode-matching,  FEM,  and  MoM 
analysis  of  microwave  circuits.  In  conclusion,  it  is 
proved  that  the  high  efficiency  and  effectiveness  of 
WBRA  turns  the  strategy  of  bandwidth  reduction 
combined  with  a  banded  solver  into  the  most  prof¬ 
itable  way  of  solving  linear  systems  in  electromag¬ 
netic  numerical  methods. 

I.  Introduction 

The  use  of  numerical  methods  is  nowadays  com¬ 
monly  accepted  as  the  most  effective  and  efficient 
way  to  attack  and  solve  electromagnetic  (EM) 
problems.  Many  numerical  codes  are  routinely  used 
in  the  CAD  of  EM  circuits,  in  complex  scatter¬ 
ing  analyses,  and  in  electromagnetic  compatibility 
evaluations,  just  to  mention  some  possible  indus¬ 
trial  and  research  tasks  which  are  currently  mainly 
performed  via  numerical  approaches. 

The  obvious  consequence  of  the  continuous 
growth  of  numerical  methods  in  daily  work  is  an 
increasing  demand  for  numerical  efficiency  and  per¬ 
formance.  As  the  problems  get  more  complex,  the 
issue  of  optimum  memory  exploitation  and  CPU¬ 
time  reduction  is  crucial,  provided  that  suitable  nu¬ 
merical  accuracy  be  guaranteed. 


A  common  bottle-neck  limiting  the  performance 
of  many  numerical  approaches  is  represented  by  the 
solution  of  linear  systems,  usually  sparse,  which 
is  very  often  one  of  the  strongest  numerical  tasks 
for  many  EM  methods.  Mode-matching  (MM), 
Method  of  Moments  (MoM),  both  in  its  standard 
formulation  and  when  using  wavelet  expansions, 
and  Finite  Element  Methods  (FEM)  are  just  some 
examples  demonstrating  this. 

In  previous  papers  [17],  [3]  it  was  demonstrated 
that  a  very  efficient  strategy  to  improve  the  per¬ 
formance  of  many  EM  codes  is  the  enhancement 
of  the  linear  system  solution  time  by  an  appropri¬ 
ate  transformation  of  the  system  matrix.  The  ma¬ 
trix,  generally  sparse,  is  transformed  into  a  banded 
one,  with  reduced  bandwidth,  this  paving  the  way 
for  a  very  effective  use  of  high-performing  banded 
solvers.  The  performance  of  the  algorithm  to  re¬ 
duce  the  bandwidth  of  a  sparse  matrix  is,  in  this 
perspective,  a  key-point.  In  this  paper,  a  new 
method  is  proposed  to  accomplish  this  task.  It  is 
suited  to  every  kind  of  sparse  matrix,  but  specifi¬ 
cally  tuned  to  achieve  maximum  performance  on 
typical  matrix  patterns  of  EM  problems.  It  is 
proved  to  outperform  all  the  previous  approaches, 
including  commercial  packages,  on  several  real  EM 
cases.  The  availability  of  such  an  efficient  band¬ 
width  reducer  turns  its  use,  in  conjunction  with 
a  banded  solver,  into  the  most  effective  solution 
strategy,  differently  from  before,  when  the  band¬ 
width  minimization  effort  was  not  so  profitable, 
and  the  use  of  a  sparse  solver  without  matrix  pre¬ 
processing  was  sometimes  the  winning  approach. 

The  paper  is  structured  as  follows.  In  section  1 
an  overview  of  the  problem  is  proposed.  In  sec¬ 
tion  2  the  new  algorithm  for  bandwidth  reduction 
is  proposed.  In  section  3  results  are  given.  Finally, 
conclusions  are  drawn. 

II.  Bandwidth  reduction  used  in 

CONJUNCTION  WITH  DIRECT  SOLVERS  VS. 

ITERATIVE  SPARSE  METHODS 

Let 

Ax  =  B  (1) 
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be  a  sparse  system. 

We  know  from  very  basic  matrix  algebra  that, 
considered  a  permutation  matrix  P  so  that 
P  •  PT  =  I  (where  I  is  the  identity  matrix),  the 
system 

PAPtx'  =  PB  (2) 

has  the  same  numerical  stability  as  (1),  as  its  condi- 
tion  number  is  not  changed.  Moreover,  the  solution 
x  of  (1)  is  easily  found  from  the  solution  x'  of  (2), 
as 

x  =  P-V  (3) 

Therefore,  if  the  permutation  matrix  is  appro¬ 
priately  chosen,  so  that  the  transformed  system 
matrix  A'  =  PAPt  is  banded,  with  small  band¬ 
width,  banded  solvers  can  be  used  with  very  high 
performance,  as  their  complexity  depends  quadrat- 
ically  on  the  matrix  bandwidth  [5].  Unfortunately, 
the  task  of  finding  the  transformation  P  that  gives 
the  minimum  possible  bandwidth  is  a  very  difficult 
problem;  in  [12]  it  is  shown  to  be  an  NP-hard  prob¬ 
lem. 

Many  efforts  have  been  made  until  now  in  or¬ 
der  to  solve  this  problem  both  in  an  exact  way 
providing  the  optimal  permutation,  and  in  an  ap¬ 
proximated  way,  so  that  good  permutations  that 
sensibly  reduce  the  bandwidth  can  be  found  with 
a  much  smaller  computational  effort.  A  very  thor¬ 
ough  review  of  approaches  is  proposed  in  [5],  which 
must  be  integrated  with  more  recent  contributions 

[4],  [7],  [2], 

Due  to  this  difficulty,  instead  of  looking  for  an 
optimum  transformation,  alternative  approximated 
but  faster  strategies  are  usually  followed  in  practi¬ 
cal  applications,  as  illustrated  in  the  next  section 
where  the  method  for  bandwidth  reduction  pro¬ 
posed  in  [6]  referred  as  WBRA  will  be  reviewed, 
putting  forward  its  peculiar  features  with  respect 
to  other  previous  approaches. 

Of  course,  the  development  of  an  effective  band¬ 
width  reduction  algorithm  does  not  guarantee  an 
a-priori  enhancement  of  performance  whatever  the 
EM  numerical  method  might  be.  Its  impact  must 
be  benchmarked  by  comparing  it  with  respect  to 
the  current  existing  approaches.  Among  them, 
a  usual  approach  is  the  use  of  iterative  sparse 
solvers,  which  are  generally  applied  directly  to  the 
original  matrix,  without  any  previous  preprocess¬ 
ing.  Iterative  solvers  are  not  very  efficient,  and  re¬ 
quire  large  numbers  of  iterations  especially  on  non¬ 
well-conditioned  problems.  Nonetheless,  they  are 
largely  used,  due  to  their  easy  availability  and  to 
their  strong  stability.  In  fact,  direct  sparse  solvers, 
which  could  be  in  principle  more  performing,  are 
often  prone  to  the  risk  of  dense  LU  factors  [5], 
highly  degrading  memory  and  computing-time  re¬ 
quirements,  and  are  therefore  not  considered  a  vi¬ 
able  solution. 


Therefore,  on  the  basis  of  the  previous  consider¬ 
ations,  when  presenting  results,  we  generally  pro¬ 
pose  a  comparison  between  a  strategy  using  a  band¬ 
width  reduction,  and  a  more  standard  one,  using  an 
iterative  sparse  solver. 

III.  The  WBRA  Approach  for  Bandwidth 
Reduction 

One  of  the  most  effective  classes  of  algorithms 
specifically  devoted  to  bandwidth  reduction  is  the 
one  derived  from  the  Cuthill-McKee  (CM)  method 
[8] .  The  main  idea  of  this  class  of  algorithms  is  re¬ 
lated  to  the  graph  representation  of  the  matrix  (see 
Fig.  1).  Consider  a  symmetrically  structured  ma¬ 
trix  A  of  order  n,  let  G  =  ( N ,  E)  be  the  undirected 
adjacency  graph  related  to  A,  where  each  node 
i  £  N  =  {1, . . . ,  n}  represents  the  2-th  row/ column 
of  the  matrix,  and  there  is  an  edge  (i,j)  between 
two  nodes  i  and  j  (i  ^  j)  if  and  only  if  the  ele¬ 
ment  of  matrix  A  aij  ^  0.  The  basic  idea  of  the 
computation  can  be  summarized  by  the  following 
steps: 

1)  partitioning  phase :  select  a  root  node 
r,  and  partition  N  into  subsets  called  levels 
{Lo,Li)...')Lp}  with  r  €  To,  so  that  there  are 
edges  only  between  nodes  belonging  to  the  same 
level  or  to  two  adjacent  levels;  a  partition  into  lev¬ 
els  is  called  level  structure.  In  Fig.  2,  the  level 
structure  obtained  with  root  equal  to  6  is  shown. 

2)  numbering  phase :  sort  the  nodes  by  increas¬ 
ing  level  index,  and  inside  each  level  number  them 
according  to  a  particular  criterion. 

As  the  other  algorithms  in  the  CM  class,  WBRA 
follows  this  general  scheme,  but  in  the  two  phases 
exploits  the  structure  of  the  combinatorial  opti¬ 
mization  problem  underneath. 

In  the  partitioning  phase  a  partition  into  lev¬ 
els,  whose  larger  subset  has  minimum  cardinality, 
is  sought.  In  fact  the  bandwidth  is  directly  af¬ 
fected  by  the  size  of  the  largest  subset.  However,  as 
this  problem  is  as  difficult  as  finding  the  minimum 
bandwidth,  an  approximated  algorithm  is  applied. 
In  Fig.  3  a  possible  redistribution  of  nodes  be¬ 
tween  levels  is  provided.  In  this  case  the  width  of 
the  largest  level  in  the  new  structure  is  reduced  to 
2. 

Now  consider  the  numbering  phase.  WBRA 
applies  the  numbering  to  a  set  of  “promising” 
level  structures  determined  during  the  partition¬ 
ing  phase,  that  is  to  a  set  of  level  structures  whose 
largest  level  is  small. 

Let  {To,  L\ , . . . ,  Lp}  be  the  partition  under  con¬ 
sideration  during  the  generic  iteration: 

Numbering  the  nodes,  we  assign  the  first  num¬ 
bers  (i.e.  the  first  positions  in  the  permuted  ma¬ 
trix)  to  the  nodes  in  Lo,  then  the  other  positions 
are  assigned  following  the  precedence  between  lev¬ 
els,  that  is  any  element  of  level  Lh-i  always  has  a 
smaller  number  with  respect  to  any  other  element 
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of  level  Lh,  h  =  1, . . .  ,p.  Thus  the  numbering  phase 
is  carried  out  by  a  sequence  of  steps,  one  for  each 
level. 

Determining  the  optimal  numbering  of  level  Lh 
WBRA  considers: 

•  (i)  the  numbers  assigned  to  the  elements  of 
Lh- 15  but  in  addition  it  considers  also 

•  (ii)  the  possible  effects  of  the  numbering  in 
level  Lh+ 1. 

These  two  criteria  give  rise  to  a  well  charac¬ 
terized  combinatorial  optimization  problem  whose 
particular  structure  allows  us  to  determine  the 
numbering  of  a  level  in  almost  linear  time.  The 
numbering  obtained  by  applying  the  algorithm  to 
the  level  structure  of  Fig.  3  is  shown  in  Fig  4.  A 
detailed  discussion  of  the  combinatorial  structure 
of  the  problem,  and  of  the  numbering  algorithm 
can  be  found  in  [6]. 

The  rearranged  matrix  according  to  the  permu¬ 
tation  found  by  the  algorithm  is  shown  in  Fig.  5. 

A.  The  unsymmetric  case 

The  method  proposed  until  now,  based  on  a  ma¬ 
trix  representation  through  an  adjacency  graph, 
is  devoted  to  matrices  with  symmetrical  patterns. 
On  the  contrary,  as  described  in  the  result  sec¬ 
tion,  a  method  working  also  on  matrices  with 
an  unsymmetric  zero-non  zero  structure  is  needed 
quite  often.  One  of  the  key-points  of  WBRA  is 
its  amenability  to  cope  with  this  problem,  with¬ 
out  degradating  the  performance  of  the  algorithm. 
This  is  achieved  by  symmetrizing  the  matrix  struc¬ 
ture  in  a  cumbersome  way. 

Generally,  the  matrix  pattern  could  be  sym¬ 
metrized  in  two  possible  fashions: 

i)  consider  a  matrix  A*  where 

ii)  symmetrize  the  matrix,  that  is  consider  A*  = 

*aa*. 

The  latter  approach  is  avoided  as  it  may  intro¬ 
duce  some  ill  conditioning.  In  the  former  case,  any 
symmetric  bandwidth  reduction  algorithm  can  be 
applied  to  A* ,  then  matrix  A  is  permuted  accord¬ 
ing  to  the  obtained  reordering.  The  drawback  of 
this  approach  is  evident  when  A  is  highly  unsym¬ 
metric:  many  zero  elements  are  dealt  as  they  would 
be  non  zero.  This  is  why  devising  ad  hoc  algo¬ 
rithms  that  can  take  advantage  of  the  unsymmetry 
becomes  important. 

By  contrast  to  the  symmetric  case,  in  the  band¬ 
width  minimization  of  unsymmetric  matrices,  we 
are  not  obliged  to  apply  the  same  permutation  to 
rows  and  columns.  Thanks  to  this  observation,  we 
propose  an  algorithm  divided  into  two  phases.  In 
the  former  phase  a  permutation  is  applied  to  rows 
only.  Then  the  matrix  is  symmetrized  according 
to  method  i).  In  the  latter  phase  the  bandwidth 
is  reduced  by  applying  the  same  permutation  to 


rows  and  columns.  Finally  the  matrix  is  “desym- 
metrized”  by  removing  the  nonzero  elements  intro¬ 
duced  by  the  symmetrization  step. 

We  adopted  two  methods  for  the  first  phase.  The 
former  method  applies  the  row  permutation  that 
maximizes  the  number  of  non  zero  elements  on  the 
diagonal.  This  problem  is  known  as  the  transversal 
maximization  [5].  An  alternative  method  ( symme¬ 
try  maximization )  tries  to  maximize  the  number  of 
symmetric  elements,  that  is  tries  to  permute  the 
matrix  rows  so  that  in  the  final  matrix  the  num¬ 
ber  of  elements  ^  0  having  the  corresponding 
aji  ^  0  are  as  many  as  possible. 

After  the  first  phase  (either  transversal  maxi¬ 
mization  or  symmetry  maximization),  the  WBRA 
is  applied.  The  final  performance  (as  demonstrated 
by  results)  is  quite  attractive.  The  algorithm  is 
publically  available  through  Internet  at  the  web  site 
http : //dvorak . istel . ing . unipg . it . 

IV.  Results 

Four  main  areas  of  applications  are  proposed: 
the  MM  analysis  of  rectangular  waveguide  circuits, 
the  FEM  analysis  of  planar  circuits  with  metalic 
boxes,  the  MoM  analysis  of  planar  circuits  with 
a  Mixed-Potential  Integral  Equation  (MPIE)  ap¬ 
proach,  and  a  MoM  using  wavelet  expansions.  In 
the  case  of  MM  and  MoM,  the  use  of  a  band¬ 
width  reduction  in  couple  with  a  banded  solver  is 
compared  with  a  standard  iterative  biconjugate- 
gradient  sparse  solver  with  preconditioning.  For 
FEM,  the  performance  of  a  standard  package  is 
compared  with  an  implementation  taking  advan¬ 
tage  from  bandwidth  reduction. 

A .  MM  Analysis  of  MW  Circuits 

The  MM  analysis  of  MW  circuits  is  an  efficient 
and  rigorous  method,  often  used  in  CAD  packages. 
Among  its  several  attractive  issues,  the  amenabil¬ 
ity  to  MW  circuit  optimization,  via  the  Adjoint 
Network  Method  is  one  of  the  most  interesting,  as 
well  as  a  paramount  impulse  to  improve  its  perfor¬ 
mance. 

As  already  discussed  and  demonstrated  in  [17], 
the  solution  of  a  sequence  of  linear  systems,  with 
different  right-hand-sides  and  same  matrix,  is  the 
numerical  core  of  the  approach.  The  pattern  of  the 
sparse  system  matrix  A  depends  on  the  numbering 
adopted  for  the  physical  and  electrical  ports,  as 
well  as  on  the  number  of  modes  selected  in  every 
section  of  the  circuit. 

In  this  paper,  we  compare  the  performance  of 
WBRA  with  other  methods  to  reduce  the  band¬ 
width  of  the  system  matrix,  and  show  the  corre¬ 
sponding  performance  of  the  MM  analysis.  More 
specifically,  we  compare  WBRA  with  a  proprietary 
routine  performing  the  Modified-  Reversed- Cut hill- 
McKee  (MRCM)  approach  (one  of  the  most  effi¬ 
cient  implementation  of  the  CM  method  [8]),  with 
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a  commercial  routine  for  bandwidth  reduction  pro¬ 
posed  by  MATLAB,  and  with  a  previous  novel 
technique  proposed  by  the  authors  in  [3],  called 
Tabu  Search  (TS).  In  Tab.  I  results  are  given  on 
three  matrices  coming  from  the  MM  analysis  of  a 
complex  4x4  Butler  matrix  (Fig.  6),  with  different 
spatial  resolutions.  In  Tab.  I  the  reader  can  note: 
the  system  size  Ny  the  bandwidth  of  the  original 
matrix  A  (IB),  the  final  bandwidth  achieved  by  the 
different  algorithms  and,  in  brackets,  the  time  nec¬ 
essary  to  evaluate  the  permutation  matrix.  Times 
are  given  in  seconds,  on  an  IBM  RS6000  250  T. 


N 

IB 

MRCM 

MATLAB 

TS 

WBRA 

116 

75 

6  (0.039) 

6  (0.062) 

6  (91.1) 

6  (0.051) 

245 

115 

83  (0.24) 

55  (0.36) 

25  (242) 

27  (0.229) 

503 

452 

78  (1.572) 

61  (1.97) 

64  (1000) 

51  (1.234) 

Tab.  I:  The  effectiveness  of  different  bandwidth 
reduction  methods  and  (times  in  brackets)  their 
efficiency. 


As  inferred  from  the  table,  TS  is  too  slow  to 
be  used  on  serial  platforms,  and  is  therefore  omit¬ 
ted.  The  WBRA  is  superior  to  both  MATLAB  and 
MRCM,  as  it  is  more  effective  and  its  time  perfor¬ 
mance  is  better  as  the  size  of  the  problem  increases. 
In  Tab.  II,  for  the  same  cases  of  Tab.  I,  the  solu¬ 
tion  times  are  shown  using:  1)  a  standard  iterative 
sparse  solver  (BCG)  with  the  original  A  matrix  2) 
a  banded  solver  on  the  original  A  matrix  Q3NT)  3) 
a  banded  solver  on  the  transformed  PAP^  matrix 
(P  evaluated  with  MATLAB  (BTM))  4)  a  banded 
solver  on  the  transformed  PAPT  matrix  (P  evalu¬ 
ated  with  WBRA  (BTW)). 


N 

BCG 

BNT 

BTM 

BTW 

116 

9.5 

27.1 

0.46 

0.46 

245 

101.1 

198 

36.3 

8.2 

503 

521 

longer  than  1000 

442 

296 

Tab.  II:  Different  simulation  times  for  a  MM  code 
on  a  4x4  Butler’s  matrix.  A  standard  iterative 
solution  (BCG)  is  compared  with  a  banded 
solution  (BNT),  a  banded  solution  with 
bandwidth  minimization  by  MATLAB  (BTM), 
and  a  banded  solution  with  bandwidth 
minimization  by  WBRA  (BTW). 

As  shown  in  Tab.  I,  the  WBRA  method  outper¬ 
forms  the  other  approaches.  Its  use  (as  evidenced 
by  Tab.  II  if  you  compare  BTW  with  respect  to 
BCG),  allows  a  speed-up  in  the  system  solution  of 
up  to  one  order  of  magnitude  with  respect  to  the 
use  of  a  standard  iterative  sparse  solver,  and  (as 
evidenced  by  comparing  BTW  and  BTM)  of  up  to 
4  times  with  respect  to  the  use  of  previous  band¬ 
width  reduction  methods. 


B.  FEM  Analysis  of  Boxed  Microstrip  Lines 

The  analysis  of  microstrip  lines  surrounded  by 
a  metallic  box  is  a  common  problem  for  the  MW 
and  EMC  community.  This  problem  has  been  at¬ 
tacked  by  using  a  public  domain  FEM  package 
called  EMAP  [10].  In  EMAP,  a  key-point  in  the 
analysis  of  the  circuit  is  the  repeated  solution  of 
a  linear  system.  In  FEM,  the  system  is  generally 
reduced  to  a  banded  structure.  Therefore,  EMAP 
uses  a  banded  solver,  and  is  quite  amenable  to  be 
interfaced  with  the  above  mentioned  modules  for 
bandwidth  reduction. 

Several  tests  have  been  performed  on  circuits 
such  as  the  one  in  Fig.  7,  for  different  substrates 
and  dimensions  of  both  the  box  and  microstrip. 
Some  of  the  results  are  shown  in  two  tables,  with 
the  same  scheme  as  for  the  MM  section:  in  Tab. 
Ill  data  report  the  efficiency  of  the  different  band¬ 
width  reduction  methods,  whilst  Tab.  IV  shows 
the  effects  of  different  bandwidth  reductions  (Mat- 
lab  (BTM),  and  WBRA  (BTW))  on  the  system  so¬ 
lution  time  using  the  standard  banded  solver  used 
in  EMAP  (BNT). 


N 

IB 

MRCM 

MATLAB 

WBRA 

291 

107 

106  (0.32) 

106  (0.48) 

62  (0.29) 

615 

156 

165  (3.75) 

163  (8.9) 

88  (1.56) 

1180 

297 

245  (24.8) 

237  (142) 

184  (5.9) 

Tab.  Ill:  A  comparison  of  performance  for 
different  bandwidth  reduction  methods  for  a  FEM 
_ linear  system. _ 


N 

BNT 

BTM 

BTW 

291 

34.1 

33.8 

8.8 

615 

243 

289 

59.6 

1180 

2048 

1621 

973 

Tab.  IV:  Solution  times  for  the  FEM  problem 
using  a  standard  EMAP  banded  solver  (BNT), 
with  respect  to  1)  standard  EMAP  solver  and 
Matlab  bandwidth  reduction  (BTM)  and  2) 
standard  EMAP  solver  and  WBRA  bandwidth 
reduction  (BTW). 

It  can  be  noted  from  Tabs.  Ill  and  IV  that  the 
superior  performance  of  WBRA,  both  for  effective¬ 
ness,  and  for  computing  times,  allows  a  speed-up 
in  the  system  solution  time  of  up  to  4  times  with 
respect  to  standard  EMAP  code. 

C.  MoM  Analysis  of  Microstrip  Circuits 

Recent  enhancements  in  the  analysis  of  planar 
circuits  with  an  MPIE  approach  using  a  closed- 
form  spatial-domain  Green’s  function  [11]  allow  a 
very  efficient  implementation  of  CAD  tools.  The 
MPIE  can  be  discretized  by  using  the  MoM,  thus 
generating  a  linear  system  whose  solution  allows 
the  evaluation  of  the  scattering  parameters  of  the 
circuit.  The  system  is  generally  dense,  but  very 
recently  it  has  been  demonstrated  that  it  can  be 
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reduced  to  a  sparse  one,  without  affecting  the  ac¬ 
curacy  of  the  simulation  [1] .  Therefore,  also  in  this 
case  the  core  of  the  numerical  effort  is  a  sparse  lin¬ 
ear  system,  and  a  bandwidth  reduction  is  worth  to 
be  performed. 

In  Tabs.  V  and  VI  results  are  shown  for  two 
circuits  reported  in  Fig.  8.  As  usual,  Tab.  V 
gives  results  concerning  the  effectiveness  and  effi¬ 
ciency  of  the  WBRA  with  respect  to  other  meth¬ 
ods  for  bandwidth  reduction.  Tab.  VI  compares 
the  solution  time  using  a  standard  iterative  sparse 
solver  (BCG)  with  respect  to  using  a  banded  solver, 
invoked  straightly  (BNT)  or  after  performing  a 
bandwidth  reduction  (using  MATLAB  (BTM)  or 
WBRA  (BTW)). 


N 

IB 

MRCM 

MATLAB 

WBRA 

220 

208 

120  (0.34) 

119  (0.51) 

72  (0.3) 

401 

310 

141  (1.23) 

136  (1.5) 

82  (0.98) 

Tab.  V:  A  comparison  of  performance  for  different 
bandwidth  reduction  methods  for  a  MPIE/MoM 
linear  system. 


N 

BCG 

BTM 

BTW 

220 

71 

14.8 

6.4 

401 

214 

49.6 

18.8 

problem  using  standard  sparse  iterative  solver 
(BCG),  with  respect  to  1)  banded  solver  and 
Matlab  bandwidth  reduction  (BTM)  and  2) 
banded  solver  and  WBRA  bandwidth  reduction 
(BTW). 


From  Tabs.  V  and  VI  it  can  be  noted  that  the 
use  of  WBRA  enhances  the  efficiency  of  the  system 
solution  by  up  to  12  times  with  respect  to  a  stan¬ 
dard  sparse  iterative  solver,  and  up  to  3  times  with 
respect  to  using  a  commercial  bandwidth  reducer 
before  invoking  a  banded  solver. 


D.  MoM  using  Wavelet  expansions 

In  the  past  few  years,  the  use  of  wavelet  expan¬ 
sions  in  the  solution  of  electromagnetic  problems 
has  become  more  and  more  frequent.  Wavelet  ex¬ 
pansions  have  been  introduced,  for  instance,  in  con¬ 
junction  with  the  Method  of  Moments  (MoM)  dis¬ 
cretization  of  integral  equations  [9],  [20],  in  order 
to  solve  scattering  problems  with  large-scale  scat¬ 
tered  (thus  containing  a  variety  of  length  scales 
with  respect  to  wavelength)  [15]- [18],  or  to  analize 
slot-apertures  [19],  microstrip  floating  line  struc¬ 
tures  [21],  as  well  as  to  study  2D  and  3D  dielectric 
structures  [13],  [14].  A  common  key-issue  for  all 
the  above  mentioned  applications  is  the  derivation 
of  very  sparse  and  well-conditioned  linear  systems, 
representing  the  numerical  core  of  MoM  approaches 
[1],  [13],  [14].  The  moment  matrix  sparsity  allows 


the  use  of  very  efficient  iterative  sparse  solvers,  and 
the  good  condition  number  guarantees  a  low  num¬ 
ber  of  iterations  to  converge,  with  a  consequent 
dramatic  improvement  of  performance. 

Up  to  now,  once  the  moment  matrix  has  been 
sparsified  using  wavelet  expansions,  it  has  been  as¬ 
sumed  that  iterative  solvers  are  the  best  way  to 
attack  the  linear  system  solution.  We  demonstrate 
here  that,  by  means  of  appropriate  matrix  trans¬ 
formations,  the  use  of  a  banded  direct  solver  in 
conjunction  with  WBRA  outperforms  the  iterative 
approach,  especially  when  non-symmetric  moment 
matrices  are  attained  after  split  testing  procedures 
in  presence  of  compact-support  functions  [15],  [16], 
[13],  [14].  It  must  be  put  forward  that  the  capabil¬ 
ity  of  dealing  with  non-symmetrical  cases,  without 
lost  of  effciency,  is  one  of  the  most  attractive  fea¬ 
tures  of  WBRA.  In  fact,  for  previous  bandwidth 
reduction  approaches,  the  only  way  to  face  non- 
symmetric  problems  was  represented  by  the  matrix 
pattern  symmetrization,  with  an  obvious  dramatic 
reduction  of  performance. 

We  refer,  for  the  proposed  results,  to  a  MoM  dis¬ 
cretization  of  a  Mixed-Potential  Integral-Equation 
formulation  for  the  analysis  of  planar  microstrip 
circuits,  as  described  in  [1].  The  MoM  matrices  are 
transformed  in  accordance  with  the  use  of  Battle- 
Lemarie  multiresolution  expansions,  as  described  in 
[15],  [16],  [13],  [14],  thus  attaining  non-symmetric 
matrices  when  splitting  and  truncations  are  per¬ 
formed  to  comply  with  boundary  conditions.  A 
double-layer  microstrip  waveguide  has  been  stud¬ 
ied,  with  different  basis  expansions,  and  different 
threshold  values  vt  have  been  applied  onto  the  mo¬ 
ment  matrices,  so  that  values  having  magnitude 
less  than  vt  per  cent  of  the  largest  entry  are  con¬ 
sidered  as  zeros.  Of  course,  different  approxima¬ 
tions  are  attained  on  varying  and  errors  have 
been  estimated  by  comparing  approximate  results 
with  the  correct  result  attained  without  any  thresh¬ 
olding.  Table  VII  and  VIII  present  results  from 
two  different  cases  of  analysis  of  a  double-layer  mi¬ 
crostrip  (see  Fig.  9),  using  different  numbers  of 
Battle-Lemarie  wavelet  functions.  Different  matrix 
sizes,  respectively  N=250  and  N=478,  have  been 
attained.  For  different  thresholds,  the  matrix  spar¬ 
sity  S,  the  approximation  error,  and  the  results  are 
reported,  from  two  different  strategies:  i)  a  banded 
solver  with  WBRA  (BTW),  ii)  an  iterative  sparse 
BCG  solver  (BCG)  (the  number  of  iterations  to 
converge  is  also  reported). 


Sol.  time  (s.) 

Vt 

sl 

Error 

BTW 

BCG  (n. iterations) 

2% 

99%”” 

5.4% 

0.08 

0.09  (8) 

iW 

96% 

2.4% 

0.12 

0.23  (17) 

0.5% 

91% 

0.7% 

0.21 

0.48  (31) 

Tab.  VII  :  1 

Results  \ 

:or  a  matrix  of  dimension 
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N=250.  Computing  times  for  a  banded  solver  and 
WBRA  strategy  versus  BCG  strategy  are  shown, 
for  different  threshold  values,  and  the 
corresponding  matrix  sparsity  S  and  solution  error 
due  to  thresholding  effects.  For  BCG  the  number 
of  iterations  needed  to  converge  is  shown  in  the 
brackets. 


I  . . 

Sol,  time  (s.) 

Vt 

S 

Error 

BTW 

BCG  (n. iterations) 

2% 

98% 

4.8% 

0.30 

0.1  (11) 

1% 

94% 

2.1% 

0.64 

1.4  (14) 

0.5% 

"88%“ 

0.6% 

0.80 

6.0  (28) 

Tab.  VIII  :  Results  for  a  matrix  of  dimension 
N=478.  Computing  times  for  a  banded  solver  and 
WBRA  strategy  versus  BCG  strategy  are  shown, 
for  different  threshold  values,  and  the 
corresponding  matrix  sparsity  S  and  solution  error 
due  to  thresholding  effects.  For  BCG  the  number 
of  iterations  needed  to  converge  is  shown  in  the 
brackets. 

As  apparent  from  Tab.  VII  and  VIII,  an  ap¬ 
propriate  value  for  thresholding  is  0.5%,  so  that 
the  approximation  error  is  smaller  than  1%.  In 
this  case,  for  N=250,  a  speed-up  of  nearly  2.3  is 
achieved,  when  using  the  BTW  strategy  with  re¬ 
spect  to  BCG,  whilst  for  N=478  a  speed-up  of 
nearly  7.5  is  observed. 

V.  Conclusions 

In  this  paper  a  new  method,  called  WBRA,  to 
perform  the  bandwidth  reduction  of  a  sparse  ma¬ 
trix  has  been  presented.  It  generally  works  for 
every  kind  of  sparse  matrix,  but  it  is  specifically 
tuned  to  achieve  maximum  performance  on  typ¬ 
ical  matrix  patterns  encountered  in  electromag¬ 
netic  numerical  problems,  as  well  as  to  manage  also 
with  non-symmetric  zero-non-zero-pattern  matri¬ 
ces  (this  being  crucial  when  attacking  some  wavelet 
problems). 

The  efficiency  and  effectiveness  of  the  method 
is  superior  to  all  the  other  commercial  and  pub¬ 
lic  domain  packages  available  on  the  market,  as 
demonstrated  on  matrices  generated  by  the  mode¬ 
matching,  FEM,  MoM/MPIE,  and  MoM/wavelet 
analysis  of  rectangular  waveguide  and  planar  cir¬ 
cuits. 

The  huge  advantages  of  the  WBRA’s  use  in 
the  analysis  of  MW  circuits  is  also  proved  for  the 
above  mentioned  problems.  It  is  demonstrated  that 
enhancements  of  one  order  of  magnitude  can  be 
achieved,  with  respect  to  the  use  of  classical  itera¬ 
tive  sparse  solvers,  by  using  WBRA  in  conjunction 
with  banded  solvers.  This  is  attained  thanks  to  the 
high  efficiency  of  WBRA,  which  reduces  the  band¬ 
width  reduction  times,  improves  the  effectiveness  of 
bandwidth  reduction,  and  substantially  decreases 
the  numerical  complexity  of  the  banded  solution. 
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Fig.  1.  The  input  matrix  (bandwidth=5)  and  the  relative 
adiacency  graph. 
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Fig.  5.  The  renumbered  matrix  (bandwidth=2). 


C  C 


A:  3  dB  Coupler 
B:  0  dB  Coupler 
C:  Phase  shifter 

Fig.  6.  The  4x4  Butler  matrix  simulated  with  the  MM 
approach. 


Fig.  7.  The  boxed  microstrip  waveguide  simulated  with  the 
FEM  approach. 
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Fig.  8.  The  circuit  simulated,  with  different  values  for  wl, 
w2,  w3,  w4  and  w5,  with  the  MoM  approach. 


Fig.  9.  The  circuit  simulated,  with  different  wavelet  basis 
functions,  with  the  MoM  approach. 
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PARALLEL  IMPLEMENTATION  OF  GALERKIN  TECHNIQUE  IN 
LARGE-SCALE  ELECTROMAGNETIC  PROBLEMS 

Dimitra  I.  Kaklamani,  Konstantina  S.  Nikita  and  Andy  Marsh 
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ABSTRACT.  An  integral  equation  formulation  in 
conjunction  with  a  parallelised  Galerkin  technique  is 
employed  to  solve  large-scale  electromagnetic  (EM) 
problems.  The  proposed  technique  is  applicable  to 
EM  structures  consisting  of  similar  conducting  or 
dielectric  parts ,  defined  as  "elements".  Coupled  inte¬ 
gral  equations  are  derived  in  the  frequency  domain , 
written  in  terms  of  the  conductivity  currents  or  the 
electric  fields  developed  on  the  conducting  or  dielec¬ 
tric  " elements "  surfaces,  respectively .  The  system  of 
integral  equations  is  numerically  solved  via  the  par¬ 
allel  computed  Galerkin  technique,  with  convenient 
entire  domain  basis  functions.  Even  for  electrically 
large  structures,  the  use  of  entire  domain  basis  func¬ 
tions  leads  to  relatively  small  order  linear  systems 
and  the  main  computational  cost  refers  to  the  matrix 
fill  rather  than  the  matrix  solution.  The  parallelisation 
introduced  to  the  computation  of  the  matrix  elements 
overcomes  the  limitation  of  using  Method  of  Moments 
at  lower  and  resonant  frequencies.  The  inherent  par¬ 
allelism  of  the  introduced  technique  allows  for  the 
results  to  be  obtained  with  minimal  additional  to  the 
sequential  code  programming  effort.  Two  indicative 
electromagnetic  compatibility  applications  are  pre¬ 
sented,  concerning  the  coupling  of  incident  waves 
with  multiple  conducting  rectangular  plates  and  the 
coupling  phenomena  occurring  in  a  multi-element 
waveguide  array  looking  into  a  layered  lossy  cylin¬ 
der.  Numerical  results  are  presented,  while  the  appli¬ 
cability/suitability  of  diverse  High  Performance  Com¬ 
puting  platforms  is  judged,  based  on  both  perform¬ 
ance  obtained  and  ease  of  code  portation. 

Keywords:  Method  of  Moments,  Code  Parallelisation, 
Electrically  Large  Structures. 

1.  INTRODUCTION 

The  analysis  and  design  of  complex  realistic  electro¬ 
magnetic  (EM)  structures  is  limited  by  the  relatively 
restricted  computational  power  of  the  conventional 
sequential  computers  -even  those  of  "main  frame"  or 
"vector  processing"  type-  as  well  as  by  the  vector  na¬ 
ture  of  die  EM  radiation  and  the  extended  power  ra¬ 
diation  in  the  environment.  The  introduction  of  mas¬ 
sively  parallel  computer  architectures  has  opened  new 
research  horizons  in  this  area  (Davidson  1990).  In¬ 
deed,  the  major  advantage  of  applying  High  Perform¬ 


ance  Computing  (HPC)  to  EM  problems  (Calalo  et  al. 
1989,  Davidson  1993,  Fijany  et  al.  1995,  Cwik  et  al. 
1997,  Lu  et  al.  1997)  is  the  reduction  of  execution 
times  of  a  given  size  of  problem  from  days/hours  to 
minutes/seconds,  enabling  the  investigation  of  prob¬ 
lems,  that  were  so  computationally  expense,  that  they 
were  practically  "unsolvable".  Pioneering  research 
work  in  such  areas  becomes  an  arduous  tedious  en¬ 
deavour.  To  fully  exploit  the  computing  power  of¬ 
fered  by  available  parallel  platforms,  the  existing  al¬ 
gorithms  based  on  diverse  numerical  techniques  must 
be  re-examined  with  emphasis  on  their  efficiency  for 
parallel  implementation. 

Focusing  on  Method  of  Moments  (MoM)  algorithms 
(Harrington  1983),  their  use  in  solving  large-scale  EM 
problems  is  mainly  restricted  by  the  extensive  com¬ 
putational  cost  in  calculating  the  kernel  elements  and 
solving  the  resulting  matrix  equation.  Aiming  at  re¬ 
ducing  the  storage  requirements  and  speeding-up  the 
solution  algorithms  on  either  von  Neumann  or  parallel 
computers,  a  variety  of  basis  functions,  discretisation 
schemes  and  solvers  have  been  employed  (Davidson 
1993,  Aksun  and  Mittra  1993,  Alanen  1991,  Coen  et 
al.  1994,  Bomholdt  and  Medgyesi-Mitschang  1988). 
Subdomain,  entire  domain  and  mixed  domain  or  hy¬ 
brid  Galerkin  expansions  have  been  used  in  the  lit¬ 
erature  (Aksun  and  Mittra  1993,  Alanen  1991,  Coen 
et  al.  1994,  Bomholdt  and  Medgyesi-Mitschang 
1988).  Subdomain  basis  functions  have  been  favored, 
due  to  their  geometric  flexibility  and  ability  to  handle 
localised  surface  features,  apertures  or  feed-point  dis¬ 
tributions.  Although  the  arising  multiple  integrals  of 
the  kernel  matrix  are  relatively  easily  evaluated,  the 
kernel  matrix  becomes  of  very  large  order  for  prob¬ 
lems  even  slightly  outside  the  resonance  region,  since 
at  least  ten  basis  functions  are  approximately  required 
per  wavelength.  Therefore,  when  employing  subdo¬ 
main  Galerkin  technique,  the  main  computational  cost 
refers  to  the  matrix  solution  and  the  choice  of  the 
proper  discretisation  scheme  and  solver  is  of  major 
importance.  Alternatively,  smooth  entire  domain  basis 
functions  can  be  employed  and,  when  successfully 
selected  for  a  specific  geometry,  can  lead  to  relatively 
small  matrix  dimensions  even  for  electrically  large 
structures.  In  this  case,  the  main  computational  cost 
refers  to  the  matrix  fill  and  the  efficient  computation 
of  its  elements  is  very  crucial. 
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In  this  context,  a  drastic  reduction  of  the  computation 
time  is  achieved  in  (Park  and  Ballanis  1997),  by  in¬ 
troducing  an  analytical  technique  to  evaluate  the  as¬ 
ymptotic  part  of  the  kernel  elements.  Furthermore, 
taking  advantage  of  the  fact  that  any  EM  scattering 
solution  can  be  analysed  into  high  and  low  frequency 
components,  a  multilevel  formulation  of  MoM  has 
been  presented  in  (Kalbasi  and  Demarest  1993),  while 
in  (Kaklamani  and  Marsh  1995),  a  parallel  computed 
MoM  technique  is  used  to  analyse  fundamental  elec¬ 
trically  large  planar  conducting  structures,  both  al¬ 
lowing  substantial  computational  savings  and  there¬ 
fore  use  of  MoM  to  higher  frequencies. 

The  present  work  also  deals  with  the  parallel  imple¬ 
mentation  of  Galerkin  technique,  providing  an  effec¬ 
tive  and  accurate  near-field  solution  to  a  system  of 
frequency  domain  coupled  integral  equations  model¬ 
ling  a  class  of  canonical  large-scale  EM  problems. 
Entire  domain  basis  functions  are  used  and  emphasis 
is  given  on  the  efficient  fill  of  the  derived  matrix.  To 
this  end,  a  parallel  algorithm  is  developed,  imple¬ 
menting  a  12-point  Gauss  quadrature  integration  rule 
to  compute  the  integrals  of  the  kernel  matrix,  arising 
from  the  Galerkin  technique  procedure.  In  section  2, 
the  parallelised  Galerkin  technique  is  presented  and 
its  applicability  to  different  EM  engineering  areas  is 
demonstrated.  Two  specific  problems  are  solved,  con¬ 
cerning  the  coupling  of  incident  waves  with  multiple 
conducting  rectangular  plates  and  the  coupling  phe¬ 
nomena  occurring  in  a  multi-element  waveguide  array 
looking  into  a  layered  lossy  cylinder.  The  former  ap¬ 
plication  also  serves  as  a  pilot  problem,  in  order  to 
demonstrate  the  inherent  parallelism  of  the  method 
that  can  be  exploited.  To  this  end,  in  section  3,  differ¬ 
ent  HPC  platforms  are  utilised  and  their  advantages 
are  determined  in  terms  of  both  performance  obtained 
and  ease  of  code  portation.  Finally,  section  4  provides 
some  concluding  remarks. 

2.  PARALLEL  COMPUTED  ENTIRE 
DOMAIN  GALERKIN  TECHNIQUE 

The  technique  presented  in  this  section  can  be  used  to 
analyse  a  class  of  electrically  large  canonical  struc¬ 
tures  consisting  of  similar  dielectric  or  conducting 
parts,  defined  as  "elements".  In  order  to  demonstrate 
its  applicability,  two  specific  electromagnetic  com¬ 
patibility  (EMC)  problems  are  treated. 

The  first  problem  deals  with  EM  scattering  from 
electrically  large  planar  conductors.  This  is  a  subject 
of  much  research,  since  plates  can  be  considered  as 
building  structures  of  more  complex  configurations 
(Alanen  1991,  Peters  and  Volakis  1988,  Kaklamani 
and  Uzunuglu  1994,  Ufimtsev  1996).  Furthermore, 
employing  Babinefs  principle,  the  complementary 
problem  of  EM  penetration  through  apertures  "cut"  on 
infinite  planar  conducting  screens  (Butler  et  al.  1978, 
Luebbers  and  Penney  1994)  can  also  be  considered, 
enlightening  many  EMC  problems.  The  specific  ex¬ 
amined  geometry  is  given  in  Fig.  1,  where  0-number 
of  perfectly  conducting  infinitesimal  thickness  rectan¬ 


gular  plates,  placed  with  identical  orientation  at  arbi¬ 
trary  positions  on  planes  parallel  to  the  xy-plane,  are 
presented. 


Fig .  1:  Q-number  of  perfectly  conducting  infinitesi¬ 
mal  thickness  parallel  rectangular  plates . 

The  second  problem  deals  with  the  analysis  of  cou¬ 
pling  phenomena  occurring  in  a  concentric  multi¬ 
waveguide  array  used  to  provide  focusing  inside  a 
layered  biological  tissue  model.  Focusing  of  EM  en¬ 
ergy  inside  biological  tissues  is  an  important  topic  in 
many  biomedical  applications  (Weiyian  and 
Shoroung  1992,  Arcangeli  et  al.  1984).  For  this  pur¬ 
pose,  until  now,  mainly  the  low  microwave  spectrum 
(100-1000  MHz)  has  been  employed  and  continuous 
wave  concepts  have  been  applied,  with  limited  suc¬ 
cess  (Chen  and  Ghandhi  1992,  Boag  et  al.  1993), 
while  recently  the  use  of  a  concentric  multi-element 
waveguide  array  and  pulsed  signals  of  short  pulse 
width  with  a  high  frequency  (9.5  GHz)  carrier  has 
also  been  reported  (Nikita  and  Uzunoglu  1996).  How¬ 
ever,  in  multi-element  arrays  significant  interaction 
exists  between  system  elements,  resulting  in  non- 
predictable  behaviour  of  this  type  of  systems,  as 
shown  in  (Nikita  and  Uzunuglu  1996)  for  concentric 
waveguide  systems  operating  at  low  microwave  fre¬ 
quencies.  In  the  present  paper,  the  investigation  of 
coupling  phenomena  occurring  in  concentric 
waveguide  arrays  operating  at  higher  frequencies  is 
enabled  by  applying  the  parallelised  Galerkin  tech¬ 
nique.  The  examined  geometry  consists  of  a  three- 
layer  cylindrical  biological  tissue  model,  irradiated  by 
Q  rectangular  aperture  waveguide  applicators  (Fig.  2). 
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Fig.  2:  Multi-waveguide  array  radiating  into  a 
three-layer  biological  tissue  model 


In  order  to  analyse  the  coupling  of  time  harmonic  EM 
fields  with  the  structures  shown  in  Figs.  1  and  2,  a 
system  of  coupled  integral  equations  is  derived  in 
terms  of  the  unknown  distribution  on  each  "element" 
(conducting  plate  or  waveguide  aperture)  surface,  by 
using  a  Green’s  function  approach  or  a  phase  space 
eigenwaves  description  approach.  The  operator  form 
of  the  derived  coupled  system  of  integral  equations  is 
converted  into  a  matrix  equation  by  a  Galerkin  tech¬ 
nique.  Since  structures  of  canonical  shape  are  consid¬ 
ered,  entire  domain  basis  functions  are  favored  (see 
detailed  formulations  and  analyses  in  (Kaklamani  and 
Marsh  1995)  and  (Nikita  and  Uzunoglu  1996)).  The 
order  of  the  resulting  kernel  matrices  depends  upon 
the  number  of  dielectric  or  conducting  "elements"  and 
mainly  upon  the  frequency,  increasing  fast  for  prob¬ 
lems  even  slightly  outside  the  resonance  region,  since 
more  basis  functions  are  required  to  accurately  de¬ 
scribe  the  unknown  tangential  electric  or  magnetic 
field  distribution  on  each  "element"  surface.  It  is  im¬ 
portant  to  emphasise  that  the  "clever"  choice  of  the 
entire  domain  basis  functions  results  into  relatively 
small  order  matrices  and  the  main  computational  cost 
refers  to  the  evaluation  of  the  kernel  matrix  elements 
rather  than  the  matrix  solution.  Working  in  the  spec¬ 
tral  domain,  the  arising  integrals  over  the  "elements" 
surfaces  are  analytically  computed,  while  the  infinite 
phase  space  integrals,  associated  to  the  Green’s  func¬ 
tion  Fourier  transformation  or  to  the  eigenwave  field 
description,  are  numerically  evaluated,  using  a  12- 
point  Gauss  quadrature  procedure  (Abramowitz  and 
Stegun  1970).  The  exact  expressions  of  the  arising 
integrals  are  given  in  the  Appendix  for  both  problems 
under  investigation.  Their  computation  constitutes  the 
most  time  consuming  part  of  the  corresponding  de¬ 
veloped  codes  and  is  easily  parallelised,  by  distribut¬ 
ing  the  integration  over  12  processors,  without  the 
need  of  inter-processor  communication.  If  the  inte¬ 
gration  path  is  divided  into  M  independent  Gauss  cal¬ 
culations,  it  can  be  seen  that,  with  minimal  program¬ 


ming  effort,  the  method  has  an  inherent  12xA/-fold 
parallelism.  It  is,  therefore,  envisaged  that  the  com¬ 
putation  times  can  be  further  reduced  by  a  factor  of 
M,  if  12xM  processors  are  available.  These  charac¬ 
teristics  allow  for  the  algorithms  to  be  ported  to  both 
shared  and  distributed  memory  machines  without  ex¬ 
tensive  programming  effort.  Nevertheless,  since  the 
contribution  to  the  final  value  of  the  computed  inte¬ 
grals  is  non-uniform  along  the  integration  path,  spe¬ 
cial  care  has  to  be  taken  in  distributing  the  corre¬ 
sponding  calculations  over  the  12x^  processors,  in 
order  to  avoid  load  unbalancing  problems.  A  detailed 
examination  of  varying  HPC  platforms  suitability  is 
presented  in  section  3. 

Indicative  EM  results  obtained  by  applying  the  pro¬ 
posed  parallel  computed  entire  domain  Galerkin  tech¬ 
nique  are  given  in  Figs.  3  and  4.  The  convergence  and 
stability  of  the  obtained  solutions  have  been  checked 
by  increasing  the  number  of  basis  functions  used  to 
describe  the  currents  on  the  plates  and  the  electric 
fields  on  the  apertures.  Considering  the  fact  that  the 
EM  field  expressions  satisfy  both  the  Maxwell’s 
equations  and  the  relevant  boundary  conditions,  it  is 
concluded  that  the  derived  results  are  self-consistent 
and  accurate  within  the  framework  of  the  approximate 
solution  of  the  system  of  coupled  integral  equations. 
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Fig.  3:  Conductivity  currents  distribution  induced 
on  a  22x22  conducting  plate  lying  in  the  xy-plane 
and  excited  by  a  y-oriented  Hertzian  dipole  placed  at 
a  distance  d  above  its  right-back  comer. 

As  far  as  the  problem  of  EM  scattering  from  electri¬ 
cally  large  conducting  plates  is  concerned  (see  Fig.  1), 
the  convergence  rate  of  the  conductivity  currents  in¬ 
duced  precisely  on  the  plates  surfaces  is  the  most  im¬ 
portant  result  concerning  the  accuracy  of  the  proposed 
method,  since  it  refers  to  near  field  quantities.  De¬ 
tailed  convergence  tests  and  comparison  with  pub¬ 
lished  data  concerning  the  single-plate  problem  are 
presented  in  (Kaklamani  and  Marsh  1995).  Focusing 
on  edge  behaviour  phenomena,  in  Fig.  3,  it  is  shown 
how  moving  a  y-oriented  Hertzian  dipole  along  the  z 
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axis  and  towards  the  back-right  comer  of  a  2/U2A 
rectangular  plate,  affects  the  conductivity  currents 
distribution  induced  on  the  plate  surface. 

As  far  as  the  problem  of  a  layered  biological  tissue 
model  irradiated  by  a  multi-waveguide  array  is  con¬ 
cerned  (see  Fig.  2),  the  strength  of  coupling  phenom¬ 
ena  can  be  analysed,  by  exciting  one  element,  and 
computing  the  TE10  mode  coefficients  coupled  to  the 

other  elements  (mutual  coupling  coefficients,  Sc )  and 
the  coefficient  of  the  reflected  TE10  mode  on  the  same 
aperture  (self  reflection  coefficient,  Sr).  The  paral¬ 
lelised  Galerkin  technique  has  been  used  to  compute 
the  coupling  parameters  in  a  30-element  TE10 
waveguide  (2x1  cm2  aperture  size)  array  placed  sym¬ 
metrically  at  the  periphery  of  a  cylindrical  tissue 
model,  16  cm  in  diameter.  The  computations  are  car¬ 
ried  out  at  the  operation  range  (1.1  ^  -  1.8  f09f0  being 
the  cut-off  frequency  of  the  TE10  mode)  of  the  system. 
The  tissue  model  consists  of  two  layers,  simulating 
bone  and  brain  tissues  and  it  is  surrounded  by  a  2  cm 
thick  lossless  dielectric  layer.  The  thicknesses  of  the 
bone  and  the  external  dielectric  layers  are  assumed  to 
be  p2~Pi  “  0.5  cm  and  p2-p2  =  2  cm,  respec¬ 
tively  (2 p3  =20  cm)  (see  Fig.  2).  The  dielectric  con¬ 
stant  of  the  external  layer  is  taken  to  be  s3  =  2. 1 .  The 

numerical  values  of  tissue  complex  permittivities  used 
in  the  calculations  are  defined  at  the  frequency  range 
of  interest  by  using  the  data  compiled  from  the  rele¬ 
vant  literature  (Gabriel  et  al.  1996).  In  Table  1,  con¬ 
vergence  patterns  are  presented  at  9.5  GHz,  in  terms 

of  the  self  reflection  coefficient  ( Sr )  and  the  mutual 
coupling  coefficients  with  neighbouring  (5^)  and 
opposite  (SZp)  applicators,  by  increasing  the  number 

Table  1.  Convergence  of  the  self  reflection  coefficient 
Sr  and  mutual  coupling  coefficients  (S„e,  Scop)  for 

the  configuration  of  Fig. 2  at  9.5  GHz  by  increasing 
the  number  of  aperture  modes. 


Modes  appearing 
on  the  excited 
aperture 

H 

u 

m 

TE10 

0.61 

Z-14.350 

0.0091 

Z-44.9° 

0.0027 

Z- 126.6° 

te10te12 

TM12 

0.612 

Z-15.2° 

0.0123 

Z-41.70 

0.0024 

Z-130.3® 

BgSjjMjSI 

iilJMMiiM 

0.6035 

Z-14.020 

0.011 

Z-42° 

0.0025 

Z-128° 

of  modes  appearing  on  the  excited  aperture.  It  can 
easily  be  observed  that  the  subset  of  modes  (TE10, 
TE12,  TE30,  TE32,  TM12,  TM32)  appearing  on  applicator 
apertures  is  considered  to  be  sufficient  to  assure  con¬ 
vergence  and  accuracy.  Numerical  results  for  the 

strength  of  coupling  between  neighbouring  (  S*e  )  and 
opposite  (S^,)  applicators  in  the  examined  30- 


element  waveguide  array  are  presented  in  Fig.  4,  at 
the  operation  range  of  the  system.  It  is  important  to 
emphasise  that  in  the  obtained  results,  a  detailed 
three-dimensional  EM  model  is  employed,  which 
takes  into  account  the  modification  of  the  field  on 
each  waveguide  aperture  resulted  from  the  other  radi¬ 
ating  elements  of  the  array,  as  well  as  from  the  pres¬ 
ence  of  the  lossy,  layered,  dielectric  body  standing  at 
the  near  field  region. 


Fig.  4:  The  ratio  of  the  mutual  coupling  coefficient 
between  neighbouring  ( Scne)  and  opposite  ( Scop)  ap¬ 
plicators  to  the  self  reflection  coefficient  (Sr)  in  a 

symmetric  configuration  of  30  rectangular  aperture 
(2x1  cm2)  waveguide  applicators,  placed  at  the  pe¬ 
riphery  of  a  cylindrical  body  model  of  circular  cross 
section,  20  cm  in  diameter,  at  the  operation  range  of 
the  system. 

3.  IMPLEMENTATION  ON  VARIOUS  HPC 
PLATFORMS 

Due  to  the  nature  of  the  computation  of  the  phase 
space  integrals,  which  are  the  elements  of  the  system 
kernel  (see  Appendix),  the  developed  algorithm  can 
be  adapted  for  both  shared-memory,  distributed- 
memory  and  additionally  vector  processing  HPC  plat¬ 
forms.  It  has  already  been  shown  in  (Kaklamani  and 
Marsh  1995)  that  the  code  possesses  a  substantial 
amount  of  inherent  parallelism,  that  can  be  exploited 
by  a  shared-memory  architecture.  Nevertheless, 
varying  HPC  platforms  can  exploit  the  different  de¬ 
grees  of  parallelism  inherent  in  the  algorithm. 

The  first  application,  shown  in  Fig.  I,  concerning  the 
EM  illumination  of  Q  number  of  conducting  rectan¬ 
gular  plates  lying  on  parallel  planes,  also  serves  as  a 
pilot  problem,  in  order  to  answer  the  emerging  ques¬ 
tion,  which  HPC  platform  to  use  in  solving  more 
complex  EM  problems.  For  reasons  of  CPU  time 
savings,  the  case  Q=  1  and  ax-bx  (i.e.  one  square 
plate)  is  chosen  for  the  benchmarking  process.  Small, 
intermediate  and  large  size  problems  (AM 8,  N- 98  and 
N=  512  respectively,  with  N  denoting  the  order  of  the 
derived  linear  system)  are  solved  on  each  HPC  plat- 
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form  and  the  corresponding  CPU  times  are  compared. 
There  are  two  crucial  points  that  should  be  noted  here. 
Firstly,  the  three  orders  of  systems  that  have  been 
selected  (A=8,  JV=98,  N=5\2)  are  only  indicative,  in 
order  to  demonstrate  the  scalability  of  the  proposed 
method.  That  is,  by  filling  and  solving  the  same  order 
system,  other  geometries  (e.g.  larger  number  of  plates 
at  a  lower  frequency)  can  be  solved.  Secondly,  the 
characterisation  of  the  three  cases  as  small,  intermedi¬ 
ate  and  large  size  problems  is  not  referred  to  the  ab¬ 
solute  value  of  the  system  order,  but  to  the  electrical 
size  of  the  problem  that  can  be  treated.  The  system 
order  is  in  any  case  small,  due  to  the  employment  of 
entire  domain  basis  functions.  For  example,  by  the 
jV=512  system,  the  case  of  one  square  plate  of  1(U 
side  can  be  accurately  solved. 

The  HPC  platforms  currently  evaluated  are  the  Silicon 
Graphics  Power  Challenge  and  Intel  Paragon,  both 
located  at  the  NTUA,  Greece  and  the  CRAY  C90  and 
CRAY  T3D,  both  located  at  CINECA,  Italy.  The 
shared-memory  Silicon  Graphics  Power  Challenge 
has  14  TFP  Processors  (MIPS  R8000  CPU  &  MIPS 
R8010  FPU),  a  16K  data  cache,  a  16K  instruction 
cache  and  8- way  interleaved  main  memory  of  1024 
Mbytes,  providing  a  peak  performance  of  4.2  Gflops, 
while  the  distributed-memory  Intel  Paragon  XP/S  has 
48  processing  nodes  with  i860  processors,  providing  a 
peak  performance  of  48  x  75  Mflops  =  3.6  Gflops. 
The  shared-memory  CRAY  vector  processing  C-90 
has  2x1  Gflops  vector  processing  CPU’s  providing  a 
peak  performance  of  2  Gflops  and  the  massively  par¬ 
allel  distributed-memory  CRAY  T3D  has  64  proc¬ 
essing  nodes  providing  a  peak  performance  of  64  x 
150  Mflops  =  9.6  Gflops.  Each  of  the  HPC  platforms 
is  judged  on  both  performance  obtained  and  ease  of 
code  portation,  which  is  directly  related  to  the  extent 
of  extra  needed  code. 

Before  comparing  these  four  computing  platforms, 
some  generalised  assumptions  need  to  be  mad.  Firstly, 
the  platforms  are  aimed  to  be  used  only  as  computa¬ 
tional  tools.  The  main  concern  is  related  to  the  ease  of 
code  portation  and  the  ability  to  solve  the  largest 
problem  in  a  given  tolerable  period  of  time,  without 
getting  involved  with  the  pragmatics  and  detailed 
code  optimisations.  Secondly,  if  a  computational  plat¬ 
form  with  a  performance  of  1  Gflop  takes  t  hours  to 
solve  a  given  problem,  then  the  same  platform  with  a 
2  Gflops  performance  would  need  t/2  hours  to  solve 
the  same  problem.  Thirdly,  since  the  platforms  are  of 
different  sizes  and  costs,  two  performance  figures  are 
of  interest:  the  absolute  time  taken  to  solve  a  given 
problem  ("absolute  performance")  and  an  estimated 
relative  time  taken  to  solve  the  same  problem  ("rela¬ 
tive  performance").  This  relative  time  is  calculated  as 
the  time  taken  if  the  platform  had  a  1  Gflop  peak  per¬ 
formance.  These  assumptions  enable  trans- 
architectural  platforms  to  be  composed  at  an  abstract 
level,  without  getting  involved  with  complex  bench¬ 
marks  etc.,  as  discussed  in  detail  in  the  following  and 
summarised  in  Table  2. 


As  far  as  the  shared-memory  SGI  Power  Challenge 
platform  is  concerned,  the  12-point  quadrature  Gauss 
algorithm  is  parallelised,  by  subdividing  the  integra¬ 
tion  path  and  splitting  the  corresponding  calculations 
onto  12  processors.  This  is  achieved,  by  the  hand  ad¬ 
dition  of  a  single  line  of  code,  containing  the  parallel 
directive  CSDOACROSS.  For  brevity  and  ease  of 
porting  the  code,  only  12  of  the  possible  14  TFP 
processors  are  used.  The  peak  performance  used  is 
therefore  12  x  300  Mflops  =  3.6  Gflops.  The  con¬ 
sumed  CPU  time  is  40  seconds  for  N=  8,  23  minutes 
for  N- 98  and  approximately  25  hours  for  A=512, 
resulting  in  relative  performances  of  approximately 
3.6  x  40  =  144  seconds,  3.6  x  1380  =  4968  seconds 
and  3.6  x  9000  =  32400  seconds  respectively  (see 
Table  2).  The  large  size  problem  CPU  time  (25  hours) 
also  includes  a  large  overhead,  due  to  page  swapping, 
incurred  by  accessing  huge  data  arrays.  Fortunately, 
the  porting  from  standard  FORTRAN  77  to  the  SGI 
Power  FORTRAN  has  only  taken  30  minutes  of  pro¬ 
gramming  and  consultation  effort.  When  employing 
more  rigorous  optimisation  techniques,  such  as  using 
the  full  potential  of  software  pipelining  and  the  SGI 
PFA  (Parallel  FORTRAN  Accelerator),  it  is  probable 
that  these  execution  times  will  be  considerably  further 
reduced. 


Table  2 .  Performance  comparison  of  diverse  HPC 
machines. 


Machine 

SGI 

Power  Challenge 

Intel  Paragon 
XP/S 

CRAY 

C-90 

CRAY  T3D 

Peak  Performance 
(Gflops) 

3.6 

2775 

2.0 

9.0 

Absolute 

A=8 

40 

47 

255 

26 

Perform. 

AH  98 

1380 

2940 

1960 

2040 

(sec) 

N~5\2 

9000 

16920 

12240 

- 

Relative 

A=8 

144 

130 

510 

234 

Perform. 

AH  98 

4968 

8158 

3920 

18360 

(sec) 

N-5\2 

32400 

46953 

24480 

- 

Due  to  the  fact  that  the  algorithm  can  be  divided  to  M 
independent  Gauss  calculations,  as  described  in  sec¬ 
tion  2,  a  greater  number  of  processors  can  be  used. 
Therefore,  the  parallelisation  of  the  algorithm  or  the 
distributed-memory  Intel  Paragon  XP/S  consists  in 
using  37,  out  of  the  possible  48,  i860  processing 
nodes;  3  groups  of  12  slave  processing  nodes  (i.e. 
M=  3  independent  Gauss  calculations)  and  a  master 
controlling  processing  node,  leading  to  a  peak  per¬ 
formance  of  37  x  75  Mflops  =  2.775  Gflops.  Com¬ 
munication  between  processing  nodes  is  accom¬ 
plished  via  explicit  message-passing.  Each  group  of 
12  processing  nodes  concentrates  on  the  parallelisa¬ 
tion  of  a  single  12-point  quadrature  Gauss  algorithm, 
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analogous  to  the  SGI  approach.  It  was  seen,  however, 
that  this  resulted  in  load  imbalances,  because  of  the 
fact  that  all  the  kernel  integrands  given  by  eq.  (A.l) 
vary  more  rapidly  with  the  increase  of  £;  therefore, 
when  integrating  with  respect  to  the  ^-variable  more 
subdivisions  must  be  considered  for  large  k  values, 
leading  to  further  computational  cost  with  the  k  in¬ 
crease.  After  an  explicit  investigation,  the  optimum 
distribution  was  derived.  To  this  end,  the  porting  from 
standard  FORTRAN  77  to  the  Paragon  platform  has 
taken  about  2  hours  involving  about  15  modifications 
of  the  original  code.  These  modifications  involve  the 
introduction  of  the  message-passing  constructs  irecv 
and  isend  and  the  explicit  subdivision  of  the  total 
computation  over  the  36  slave  processing  elements. 
The  consumed  CPU  time  is  47  seconds  for  N=  8,  49 
minutes  for  A=98  and  approximately  47  hours  for 
A=512,  resulting  in  relative  performances  of  2.775  x 
47  «  130  seconds,  2.775  x  2940  «  8158  seconds  and 
2.775  x  16920  =  46953  seconds  respectively  (see  Ta¬ 
ble  2).  As  it  is  expected,  due  to  the  peak  performances 
offered,  the  Paragon  performs  slower  than  the  SGI 
machine.  However,  what  is  more  disappointing  is  the 
relative  performances  of  the  Paragon  compared  to  the 
SGI  machine.  For  example,  when  N=  98,  the  relative 
performance  of  the  Paragon  machine  is  81158  sec¬ 
onds  compared  with  the  SGI  machine  relative  per¬ 
formance  of  4968  seconds.  In  this  case,  for  this  ex¬ 
ample,  a  Paragon  machine  with  1.6  times  the  peak 
performance  of  the  SGI  machine  would  be  required  to 
solve  the  given  problem  in  the  same  time.  Clearly 
some  work  has  to  be  done  to  optimise  the  data  local¬ 
ity,  when  using  the  Paragon  to  solve  problems  with 
very  large  data  arrays.  It  must  be  noted  that,  when 
N= 98,  the  message-passing  and  system  overheads 
represents  about  2%  of  the  total  execution  time, 
whereas,  when  A=512,  this  value  increases  to  25%. 
From  these  initial  investigations,  it  can  be  seen  that, 
the  algorithm  possesses  an  inherent  coarse-grained 
parallelism  that  can  be  exploited  by  both  architectural 
models. 

The  suitability  of  the  CRAY  C-90  architecture  is  also 
investigated,  to  see  if  the  algorithm  possesses  a  sub¬ 
stantial  fine-grain  parallelism,  that  could  be  exploited 
by  a  vector  processor.  Each  of  the  two  CPUs  of  the 
CRAY  C-90  is  a  vector  processor  with  128  banks,  a 
clock  of  4  ns  and  a  memory  speed  of  88  ns.  The  ex¬ 
ploitation  of  vectorisation  and  parallelisation  is 
achieved  by  the  CRAY  cf77  compiler  with  default 
settings  and  with  aggressive  optimisation  switched  on. 
The  addition  of  two  vector  processing  directives 
(CDIRS  IVDEP)  leads  to  further  reduction  of  the  exe¬ 
cution  times.  The  porting  process  required  about  1 
hour.  It  consisted  mainly  in  examining  the  code  listing 
for  further  vectorisable  loops  and  introducing  2  line 
modifications.  The  consumed  CPU  time  is  255  sec¬ 
onds  for  N=  8,  33  minutes  for  N=  98  and  approxi¬ 
mately  34  hours  for  jV=512,  resulting  in  relative  per¬ 
formances  of  2.0  x  255  =  510  seconds,  2.0  x  1960  = 
3920  seconds  and  2.0  x  12240  =  24480  seconds  re¬ 


spectively  (see  Table  2).  This  indicates  that  the  algo¬ 
rithm  contains  a  significant  amount  of  fine-grain  par¬ 
allelism  that  can  be  exploited  by  the  vector  processor. 
The  relative  performances  are  impressive,  out  per¬ 
forming  the  SGI  platform.  For  example,  from  the 
relative  performances  for  N=98,  it  can  be  seen  that  the 
CRAY  C-90  with  2.0  x  1  Gflop  =  2  Gflops  peak  per¬ 
formance  would  be  equivalent  to  a  SGI  with 
4968/(3920/2.0)  =  2.53  Gflops  peak  performance  (9  x 
300  Mflops  processors).  Considering  the  ease  of 
portability  and  absolute  performance  of  the  CRAY  C- 
90,  it  remains  a  realistic  alternative  shared-memory 
platform.  However,  both  CRAY  C-90  and  SGI  plat¬ 
forms  have  the  memory  access  bottleneck,  when  con¬ 
sidering  larger  problem  sizes. 

An  alternative  to  the  Intel  Paragon  distributed- 
memory  platform  is  the  CRAY  T3D,  a  superscalar 
architecture  with  64  processing  nodes.  Each  node  is 
based  on  the  DEC  a  21064  chip  providing  a  peak 
performance  of  150  Mflops,  with  a  clock  of  6.67ns. 
The  limiting  factor,  however,  is  the  8  Kbytes  data 
cache  and  the  8  Kbytes  instruction  cache.  Although 
possessing  facilities  for  PVM  and  message-passing, 
the  work  sharing  on  shared  data  paradigm  was  used. 
To  this  end,  porting  from  standard  FORTRAN  77  to 
the  MPP  FORTRAN  has  required  about  30  minutes 
programming  effort,  involving  definition  of  a  global 
(shared)  array  to  collect  the  results  and  explicit  subdi¬ 
vision  of  the  computation,  similar  to  the  Paragon  ap¬ 
proach.  For  brevity,  60  processing  nodes  are  used, 
providing  a  peak  performance  of  60  x  150  Mflops  = 
9.0  Gflops,  divided  in  5  groups  of  12  processing 
nodes.  Each  group  of  processing  nodes  is  responsible 
for  a  Gauss  calculation,  similar  to  that  implemented 
on  the  Paragon.  Due  to  load  balancing,  one  group  of 
processors  performs  two  Gauss  calculations.  The  re¬ 
sulting  performance  is  only  26  seconds  for  N=  8  and 
34  minutes  for  N= 98,  while  the  estimated  perform¬ 
ance  for  A=512  is  10  hours  (see  Table  2).  The  CRAY 
T3D  could  be  used  for  benchmarking  purposes  only 
for  up  to  8  hours.  Therefore,  there  are  no  available 
results  for  the  absolute  performance  in  solving  the 
N-512  problem  and  the  estimated  performance  of  10 
hours  is  considered  to  be  adequate  for  the  current  re¬ 
search.  It  must  be  noted  that,  these  results  are 
achieved  without  any  optimisation  and  there  exists  a 
load  imbalance.  Taking  these  factors  into  considera¬ 
tion,  the  performances  are  impressive. 

It  can  be  concluded  that,  for  both  ease  of  portation 
and  resulting  performance  the  shared-memory  Silicon 
Graphics  (and  to  a  lesser  extent  the  CRAY  C-90)  and 
the  distributed-memory  CRAY  T3D  appear  to  be 
more  suitable  in  solving  problems  in  the  domain  of 
electrically  large  EM  structures.  However,  the  poten¬ 
tial  of  the  shared-memory  platforms  for  solving  even 
larger  problem  sizes  is  limited  by  the  memory  access, 
whereas  the  distributed-memory  model  is  modularly 
extendible  and,  in  the  authors’  opinion,  more  suitable 
in  solving  even  larger  problem  sizes.  It  must  also  be 
noted  that  the  work  sharing  paradigm  adopted  by  the 
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CRAY  T3D  provides  an  appealing  alternative  to  mes¬ 
sage-passing  approaches,  making  programming  large 
distributed-memory  machines,  to  solve  large  problem 
sizes,  a  realistic  possibility. 


(A.i) 

where  (q=l,2,. ..,Q  ;  £=l,2,...,Q), 


4.  CONCLUSIONS 

HPC  has  been  employed  in  order  to  extent  MoM  for 
the  treatment  of  a  certain  class  of  large-scale  EM 
problems.  Namely,  a  parallel  computed  Galerkin 
technique  in  conjunction  with  a  frequency-domain 
coupled  integral  equation  formulation  has  been 
adopted,  applicable  to  electrically  large  structures 
consisted  of  similar  conducting  or  dielectric  parts. 

The  resulting  algorithm  parallelisation  overcame  the 
limitation  of  using  MoM  at  lower  and  resonant  fre¬ 
quencies,  while  the  inherent  parallelism  of  the  intro¬ 
duced  technique  allowed  for  the  results  to  be  obtained 
with  minimal  additional  to  the  sequential  code  pro¬ 
gramming  effort.  Namely,  the  phase  space  integrals, 
appearing  in  the  system  kernel  have  been  computed 
numerically  employing  a  12-point  quadrature  Gauss 
algorithm,  which  has  been  parallelised,  by  subdivid¬ 
ing  the  integration  path  and  splitting  the  correspond¬ 
ing  calculations  over  various  HPC  platforms.  For  both 
ease  of  portation  and  resulting  performance  the 
shared-memory  Silicon  Graphics  and  the  distributed- 
memory  CRAY  T3D  appear  to  be  more  suitable, 
though  the  potential  of  the  shared-memory  platforms 
for  solving  even  larger  problem  sizes  is  limited  by  the 
memory  access,  whereas  the  distributed-memory 
model  is  modularly  extendible  and,  therefore,  more 
suitable  in  solving  even  larger  problem  sizes. 

Two  specific  EMC  applications  were  solved,  con¬ 
cerning  the  coupling  of  incident  waves  with  multiple 
conducting  rectangular  plates  and  the  coupling  phe¬ 
nomena  occurring  in  a  multi-element  waveguide  array 
looking  into  a  layered  lossy  cylinder,  while  numerical 
results  were  computed  for  indicative  cases.  Due  to  the 
algorithm  parallelisation,  the  computation  times  be¬ 
came  tolerable,  allowing  the  problems’  size,  hence  the 
accuracy,  to  be  increased.  The  convergence  rate  of 
near-field  quantities  excited  on  electrically  large 
structures  is  the  most  important  result  concerning  the 
accuracy  of  the  proposed  method,  while  following  the 
same  approach,  more  complex  configurations  can  be 
constructed  and  analysed. 

5.  APPENDIX 

In  solving  the  Q  number  of  conducting  rectangular 
plates  structure  shown  in  Fig.  1,  the  analytical  expres¬ 
sions  of  the  system  kernel  sub-matrices  are  defined  as 


Kiq 

— nmn'm’ 


mo 

8**0  ak 


■At o  sgn(z°-z°Xz?-z°+<5()] 


<?0 


-j[go  *?n(-kx(x1-x0q)-ky(yl ->■“)) 
?0 


+00  +00  +C0  27r 

jj< ik_=  \dkx  \dky  •  fdkk  fd<pk  , 

Clk  -00  -00  0  0 

k_  =  kxx+kyy  =  k{cos(pkx+sw.<pky) ,  (A.2) 


kQ  =  eo^fs^/T0  is  the  free  space  propagation  constant, 
co  is  the  angular  frequency,  s0  and  //0  are  the  free 
space  dielectric  permittivity  and  magnetic  permeabil¬ 
ity  respectively,  r°g  -x°x+y^y  +  z°z  denotes  the 

position  vector  of  the  9-th  plate  centre  of  gravity  with 
respect  to  a  global  system  of  coordinates  (x,y,z), 

8 1  0+ , 

9o  =  9o(^)  = 

with  Re(#0)X)  &  lm(<yc)<0,  (A.3) 

as  required  by  the  expijcof)  time  dependence  and  the 
satisfaction  of  the  radiation  condition  at  infinity. 


A( k_)  =  A(k,cpk)  =  (kl~  k\  )xx  -  kxky(  xy  +  yx)  + 


(kl-K)yy, 


(A.4) 
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- — yy,  (A.5) 


with  Jv(x)  denoting  the  Bessel  function  of  the  first 

kind  of  v-th  order  and  S=agxbg,  (q=l,2,...,Q)  denoting 
the  finite  zero-thickness  q-ih.  conducting  plate  surface. 
The  kernel  elements  encountered  in  the  problem  of 
Fig.  2  are  of  the  following  type: 

K%  =  \fadyjfrdy'hp>t  (x,  y)  y /*',/)  • ent  (x',  /) 

s,  s, 

where  Se  is  the  aperture  of  the  ^-th  waveguide, 
en  t  and  hn  t  are  the  transverse  electric  and  magnetic 
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fields  of  the  waveguide  normal  modes,  respectively, 
being  of  TE  or  TM  type  and 


wit*1  z^(aip)=Jm{aip)  and  Z^(aip)=Ym(aip), 
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with  yn ,  Xn  being  the  propagation  constants  of  the 
TE  and  TM  modal  fields,  respectively  and 
k3  =  k0  ~e ^ .  The  matrices  involved  in  (A.7)  are 
given  by  the  following  equations: 


L  (m, t,p3)= T^l (k- p3 ) + ■ -v £1  fe fh )g; 


The  matrices  involved  in  (A.  10)  are 
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The  matrices  Z)/m  ,7V  ,i  =  l,2  and  =  1,2 ,  appear¬ 
ing  in  (A.ll),  are  obtained  from  eq.  (A. 8)  and  eq. 
(A.9),  respectively,  for  /=1,2. 
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CALL  FOR  PAPERS 


THE  APPLIED  COMPUTATIONAL 
ELECTROMAGNETICS  SOCIETY 

ANNOUNCES  A  SPECIAL  ISSUE  OF  THE  ACES  JOURNAL  ON 

EVOLUTIONARY  METHODS  IN  COMPUTATIONAL 
ELECTROMAGNETICS 

The  Applied  Computational  Electromagnetics  Society  is  pleased  to  announce  the  publication  of  a 
2000  Special  Issue  of  the  ACES  Journal  on  applications  and  advances  in  evolutionary  computing 
methods  applied  to  computational  electromagnetics.  The  objectives  of  this  special  issue  are  to 
present  applications  and  advances  in  evolutionary  methods  applied  to  antennas,  scattering,  EMC, 
microwave  filters,  and  other  relevant  electromagnetic  problems.  Prospective  authors  are 
encouraged  to  submit  papers  of  archival  value  that  address  these  objectives  and  other  suggested 
topics  listed  below. 

SUGGESTED  TOPICS 

•  genetic  algorithms 

•  genetic  programming 

•  evolutionary  algorithms 

•  electromagnetic  applications  of  computational  evolutionary  methods 

•  advances  in  GA  operators  and  algorithms 

•  advances  in  EA  and  genetic  programming 

DEADLINE  FOR  PAPERS  IS  DECEMBER  20, 1999 

Potential  contributors  wishing  to  discuss  the  suitability  of  their  contribution  to  the  special  issue 
may  contact  one  of  the  following  two  Guest  Editors  by  email  or  phone: 

Prof.  Randy  L.  Haupt  haupt@ieee.org  Tel:  (435)  797-2841 
Prof.  J.  Michael  Johnson  jmiohnson@ieee.org  Tel:  (775)  784-6485 

All  submissions  for  this  special  issue  should  be  addressed  to: 

Special  Issue  on  Evolutionary  Methods 
Prof.  Randy  Haupt 
Utah  State  University 

Dept,  of  Electrical  and  Computer  Engineering 
4120  Old  Main  Hill 
Logan,  UT  84322-4120 


CALL  FOR  PAPERS 


THE  APPLIED  COMPUTATIONAL 
ELECTROMAGNETICS  SOCIETY 

ANNOUNCES  A  SPECIAL  ISSUE  OF  THE  ACES  JOURNAL  ON 

COMPUTATIONAL  ELECTROMAGNETIC  TECHNIQUES 
IN  MOBILE  WIRELESS  COMMUNICATIONS 


The  Applied  Computational  Electromagnetics  Society  is  pleased  to  announce  the  publication  of  a 
Special  Issue  of  the  ACES  Journal  on  the  role  that  computational  electromagnetics  plays  in  the 
design  and  analysis  of  wireless  communications  components,  subsystems,  and  systems.  The 
complexity  of  present  wireless  mobile  architectures  and  the  future  complexities  of  wideband 
wireless  mobile  systems  has  sparked  an  interest  in  computational  electromagnetics  (CEM)  as  one 
of  the  many  tools  needed  for  the  design  of  third  generation  mobile  wireless  communications 
systems.  The  purpose  of  this  special  issue  is  to  draw  analysts,  designers,  and  management  from 
both  industry  and  academia  to  outline  their  ideas  and  research  on  the  role  of  CEM  in  present  or 
future  wireless  designs.  Applications  oriented  papers  are  highly  encouraged. 


SUGGESTED  TOPICS 


Applications  of  computational  electromagnetic  techniques  on  any  of  the  following: 


Smart  and  adaptive  antennas 
PCS,  Mobile,  Gateways,  Satellite  antennas 
Propagation  Models 
Atmospheric  Models 
Systems  design 

Correlation  of  measurement  techniques 


•  Electromagnetic  Interference 

•  Bioelectromagnetics 

•  LEO,  MEO,  Satellites  Communications 

•  Digital/Analog  components  design 

•  RF  components  design 
and  models 


DEADLINE  FOR  PAPERS  IS  MARCH  28,  2000 

Expected  Publication  Date  is  Fall  2000  Issue  of  ACES  Journal 

Please  submit  4  copies  of  papers  to  either  of  the  Guest  Editors  listed  below.  The  review  process 
will  commence  as  papers  are  received.  Notification  of  accepted  papers  will  be  made 
immediately  as  papers  are  reviewed. 


Ray  Perez 

c/o:Lockheed  Martin  Astronautics 
MS:  S8800  P.OBox  179 
Denver,  Colorado  80201,  USA 
phone:  303-977-5845 
fax:  303-971-4306 
emailrray.j  .perez@lmco.com 


Chris  Holloway 

NTIA/ITS.T 
325  Broadway 

Boulder,  Colorado  80303,  USA 
phone:  303-497-6184 
fax:  303-497-3680 
email:cholloway@its.bldrdoc.gov 
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THE  APPLIED  COMPUTATIONAL 
ELECTROMAGNETICS  SOCIETY 

ANNOUNCES  A  SPECIAL  ISSUE  OF  THE  ACES  JOURNAL  ON 

COMPUTATIONAL  BIOELECTROMAGNETICS 

The  Applied  Computational  Electromagnetics  Society  is  pleased  to  announce  the  publication  of  a 
Special  Issue  of  the  ACES  Journal  on  applications  and  advances  in  methods  and  applications  in 
computational  bioelectromagnetics.  The  objectives  of  this  special  issue  are  to  present  advances 
in  computational  techniques,  reviews  and/or  comparisons  of  methods,  and  applications  of 
computational  bioelectromagnetics.  Prospective  authors  are  encouraged  to  submit  papers  of 
archival  value  that  address  these  objectives  and  other  suggested  topics  listed  below. 

SUGGESTED  TOPICS 

•  Applications  of  computational  bioelectromagnetics 

•  Cellular  telephone  analysis,  design,  etc. 

•  Medical  imaging 

•  EM  Safety  analysis 

•  Etc. 

•  Methods  used  for  computational  bioelectromagnetics 

•  Finite-difference  time-domain 

•  Finite  element  method 

•  Other  methods 

•  Models  for  computational  bioelectromagnetics 

•  High  resolution  human  body  models 

•  Electrical  properties  of  human  tissue 

•  Comparisons  of  methods,  models,  or  techniques 

DEADLINE  FOR  PAPERS  IS  AUGUST  25, 2000 

Potential  contributors  wishing  to  discuss  the  suitability  of  their  contribution  to  the  special  issue 
may  contact  one  of  the  following  three  Guest  Editors  by  email  or  phone: 

Cynthia  Furse  fu  rse@ece.u  su.edu  Tel:  (435)  797-2870 

Susan  Hagness  hagness@engr.wisc.edu  Tel:  (608)  265-5739 

Ulrich  Jakobus  jakobus@ihf.uni-stuttgart.de  Tel:  +49  (0)711  685-7420 

All  submissions  for  this  special  issue  should  be  addressed  to: 

Cynthia  Furse  —  Special  Issue  on  Bioelectromagnetics 
Dept,  of  Electrical  and  Computer  Engineering 
Utah  State  University 
Logan,  UT  84322-4120 
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MTT  IMS’  00 

SESSIONS  ON  HF/VHF/UHF  TECHNOLOGY 

Boston,  Massachusetts  -  June  2000 


The  IEEE  Microwave  Theoty  and  Techniques  (MTT)  Society  has  formed  a  new  technical 
committee  (MTT- 17)  for  HF,  VHF,  and  UHF  Technology.  Its  purpose  is  to  help  MTT  address 
the  needs  of  the  26,000  RF  engineers  who  work  at  frequencies  below  1  GHz. 

We  have  sponsored  special  sessions  at  IMS’97,  and  ’98  and  a  regular  session  at  IMS’99.  We  are 
looking  forward  to  continuing  these  sessions  on  HF/VHF/UHF  technology  at  IMS ’00  in  Boston 
(June  11-16,  2000).  Papers  can  address  any  area  of  technology  that  is  of  particular  interest  to  RF 
engineers,  including: 

•  Receiver  and  associated  DSP 

•  Transmitters  and  power  amplifiers 

•  Applications  such  as  plasma/laser  drivers  and  MRI 

•  Filters  and  matching  networks 

•  Components  (transistors,  diodes,  mixers,  etc.) 

DEADLINE  FOR  PAPERS  IS  NOVEMBER  29, 1999  FOR  PAPER  FORM . 
DEADLINE  FOR  PAPERS  IS  DECEMBER  6, 1999  FOR  ELECTRONIC  FORM 

Proposals  including  a  30-50  word  abstract  and  a  500-1000  word  summary  must  be  submitted  by 
DECEMBER  6, 1999  in  electronic  form  or  by  NOVEMBER  29, 1999  in  paper  form.  For  full 
information,  see:http://www.mtt.org/ims2000.  When  submitting,  please  note  that  HF/VHF/UHF 
Technology  is  TOPIC  CATEGORY  8  on  the  author’s  information  form.  To  ensure  that  your 
submission  is  not  accidentally  routed  to  another  committee,  it  is  recommended  that  you  send  a 
copy  of  your  proposal  to  F.H.  Raab  at  the  address  shown  below. 

Please  distribute  this  notice  to  your  colleagues. 

Questions  about  the  HF/VHF/UHF  committee  or  sessions  should  be  directed  to: 

Dr.  Frederick  H.  Raab 
Green  Mountain  Radio  Research 
50  Vermont  Avenue 
Colchester,  VT,  05446,  USA 
Fax:  (802)  655-9670 
E-mail:  f.raab@ieee.org 


EUROEM  2000,  EDINBURGH,  SCOTLAND 

CALL  FOR  PAPERS 

The  Organising  Committee  has  great  pleasure  in  inviting  you  to  submit  papers  for  the  EUROEM 
2000  conference  being  held  in  Edinburgh,  Scotland  from  30  May  -  2  June  2000. 

EUROEM  2000  continues  the  tradition  of  the  EUROEM/AMEREM  Conference  Series,  drawing 
together  the  12th  High  Power  Electromagnetics  Conference,  the  5th  Ultra- Wideband  Short  Pulse 
Electromagnetics  Conference  and  the  5  Unexploded  Ordnance  Detection  and  Range  Remediation 
Conference. 

Papers  are  solicited  which  describe  original  work  suitable  for  the  three  conferences  comprising 
EUROEM  2000.  The  deadline  for  receipt  of  abstracts  is  Friday,  14  January  2000. 

Authors  are  requested  to  submit  a  one  page  abstract,  original  plus  3  copies  in  camera-ready  format 
by  the  deadline  date  to:  EUROEM  2000,  Concorde  Services  Ltd,  Suite  325,  Pentagon  Business 
Centre,  Washington  Street,  Glasgow  G3  8AZ,  Scotland,  UK.  You  should  specify  at  the  time  of 
submission  whether  your  paper  is  for  HPEM,  UWB-SP  or  UXO  and  specify  a  topic  (see  the  full  list 
on  the  website). 

The  abstract  must  consist  of  at  least  250  words  and  must  be  limited  to  one  page,  including  figures. 
Since  there  will  be  a  reduction  of  about  70%  of  the  original  linear  dimensions,  letters  and  symbols  in 
all  diagrams  should  be  sufficiently  large  and  clear.  Do  not  include  list  of  references;  a  few  open- 
literature  references  may  be  included  parenthetically,  for  example  (R  L  Lewis  &  J  R  Johler,  Radio  Sci 
2,  75-81,  1976).  Acknowledgement  of  financial  support  is  not  deemed  appropriate.  Please  note  that 
it  is  the  authors'  responsibility  and  not  that  of  the  Conference  Committee,  to  see  that  their 
abstract  and  paper  are  cleared  for  public  release. 

All  abstracts  must  be  written  in  English.  The  text  should  be  typed  single  space  on  A4  (8.27"  x  1 1.69" 
or  210  x  297  mm)  white  paper.  Margins  should  be  1.5"  on  left  and  right  and  sides  and  1"  top  and 
bottom  of  page.  The  title  should  be  centred  in  capital  letters  1  inch  from  the  top  of  the  page.  Author’s 
name  and  complete  organisational  affiliation  should  be  two  lines  below  the  title.  For  multiple 
authors,  the  expected  presenter  should  be  indicated  by  an  asterisk.  The  text  should  start  three  lines 
below  the  last  name.  Double  space  between  paragraphs. 

Alternatively,  you  may  submit  your  abstract  via  the  web.  NB:  Web  submissions  should  be  plain 
text  only,  no  graphics  allowed.  If  using  graphics  please  submit  camera-ready  copy  by  post. 

A  camera-ready  paper  will  be  requested  prior  to  the  conference,  for  accepted  UWB-SP  presentations. 
Instructions  will  be  sent  to  the  relevant  authors. 

For  additional  information  please  contact  the  Conference  Secretariat,  Concorde  Services  at  (Tel)  +44 
(0)141  221  5411  or  (Fax)  +44  (0)141  221  2411  or  visit  our  web  site  at: 

http://www.mcs.dundee.ac.uk:  8080/~euroem 


CALL  FOR  PAPERS 


.  :  ... 

/ '  A:;\  ‘  ■■  <  JESS#*. 

••  &***»  '5SP 

w  <u**t*oeom»vT*r<x*rt  j 


THE  APPLIED  COMPUTATIONAL 
ELECTROMAGNETICS  SOCIETY 


The  16th  Annual  Review  of  Progress  in  Applied  Computational  Electromagnetics 

March  20  -  25, 2000 

Naval  Postgraduate  School,  Monterey,  California 
Share  Your  Knowledge  and  Expertise  with  Your  Colleagues 

The  Annual  ACES  Symposium  is  a  large  gathering  of  engineers  and  scientists  sharing  information  and  experiences  about  the 
practical  application  of  computational  electromagnetics  (CEM).  The  symposium  offerings  include  technical  presentations,  demonstrations, 
vendor  booths,  short  courses,  and  hands-on  workshops.  All  aspects  of  CEM  are  represented.  This  Symposium  is  a  highly  influential  outlet 
for  promoting  awareness  of  recent  technical  contributions  to  the  advancement  of  CEM.  Attendance  and  professional  program  paper 
participation  from  non-ACES  members  and  from  outside  North  America  are  encouraged  and  welcome.  Papers  may  address  general  issues 
in  applied  computational  electromagnetics,  or  may  focus  on  specific  applications,  techniques,  codes,  or  computational  issues  of  potential 
interest  to  the  ACES  membership. 


Areas  and  topics: 

Code  performance  analysis 
Computer  hardware  issues 
Applications  and  Techniques: 

Antennas 

Visualization 

Integral  &  Differential  Equations 
Diffraction 


Computational  studies  of  physics 
Practical  code  application 

Bioelectromagnetics 
Wave  propagation 
Finite  Difference  &  Finite  Element 
Modal  Expansion 


Code  validation 

New,  fixed,  and  enhanced  codes 

Scattering 
Radar  cross  section 
Frequency-domain  &  Time-domain 
Physical  Optics 


INSTRUCTIONS  FOR  AUTHORS  AND  TIMETABLE 

Submission  Deadline  -  November  1,  1999:  Electronic  submission  preferred  (Microsoft  Word).  Otherwise  submit  three  copies  of  a  full- 
length,  camera-ready  paper  to  the  Technical  Program  Chairman.  Specific  format  required;  please  request  instructions  via  e-mail  or  see  on-lme 
instructions.  Authors  notified  of  acceptance  by  December  1, 1999. 

For  all  questions  regarding  the  ACES  Symposium  please  contact  Douglas  H.  Werner,  Technical  Program  Chair 
The  Pennsylvania  State  University,  21 1A  Electrical  Engineering  East,  University  Park,  PA  16802 
Tel:  (814)  863-2946,  Fax:  (814)  865-7065,  E-mail:  aces@engr.psu.edu 
or  visit  ACES  on  line  at:  http://aces.ee.olemiss.edu 


EARLY  REGISTRATION  FEES 

ACES  MEMBERS  $300  STUDENT/RETIRED/UNEMPLOYED  $130  (no  proceedings) 

NON-MEMBER  $350  STUDENT/RETIRED/UNEMPLOYED  $165(includes  proceedings) 


$500  BEST-PAPER  PRIZE 

A  $500  prize  will  be  awarded  to  the  authors  of  the  best  non-student  paper  accepted  for  the  16th  Annual  Review. 

$200  BEST  STUDENT  PAPER  CONTEST 

This  will  be  for  the  best  student  paper  accepted  for  the  16th  Annual  Review.  The  prize  will  include:  1)  free  Annual  Review  registration  for 
the  following  year;  2)  one  free  short  course;  and  3)  $200  cash  for  the  paper. 


2000  ACES  Symposium  Sponsored  by:  ACES,  NPS,  PSU,  UNR,  U  OF  NV,  MSU,  UWI,  SWRI 

In  cooperation  with:  IEEE  Antennas  and  Propagation  Society,  IEEE  Electromagnetic  Compatibility  Society  and  USNC/URSI 
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THE  APPLIED  COMPUTATIONAL 
ELECTROMAGNETICS  SOCIETY 

The  16th  Annual  Review  of  Progress  in  Applied  Computational  Electromagnetics 

March  20  -  25,  2000 

Naval  Postgraduate  School,  Monterey,  California 
Share  Your  Knowledge  and  Expertise  with  Your  Colleagues 

The  Annual  ACES  Symposium  is  an  ideal  opportunity  to  participate  in  a  large  gathering  of  EM 
analysis  enthusiasts.  The  purpose  of  the  Symposium  is  to  bring  analysts  together  to  share  information  and 
experience  about  the  practical  application  of  EM  analysis  using  computational  methods.  The  symposium 
offerings  include  technical  presentations,  demonstrations,  vendor  booths,  short  courses,  and  hands-on 
workshops.  All  aspects  of  electromagnetic  computational  analysis  are  represented. 

The  ACES  Symposium  is  a  highly  influential  outlet  for  promoting  awareness  of  recent  technical 
contributions  to  the  advancement  of  computational  electromagnetics.  Attendance  and  professional  program 
paper  participation  from  non- ACES  members  and  from  outside  North  America  are  encouraged  and  welcome. 

Papers  may  address  general  issues  in  applied  computational  electromagnetics,  or  may  focus  on 
specific  applications,  techniques,  codes,  or  computational  issues  of  potential  interest  to  the  Applied 
Computational  Electromagnetics  Society  membership. 

Areas  and  topics 

Computational  studies  of  basic  physics  Computer  hardware  issues 

Examples  of  practical  code  application  Code  validation 

New  codes,  algorithms,  code  enhancements,  and  code  fixes  Code  performance  analysis 

Partial  list  of  applications 

Communications  systems 
Remote  sensing  &  geophysics 
Dielectric  &  magnetic  materials 
Non-destructive  evaluation 
Propagation  through  plasmas 

Partial  list  of  techniques 

Diffraction  theories  Moment  methods  Physical  optics 

Frequency-domain  &  Time-domain  techniques  Hybrid  methods  Modal  expansions 

Finite  difference  &  finite  element  analysis  Numerical  optimization 

Integral  equation  &  differential  equation  techniques  Perturbation  methods 

INSTRUCTIONS  FOR  AUTHORS  AND  TIMETABLE 

Submission  Deadline  -  November  1, 1999:  Electronic  submission  preferred  (Microsoft  Word).  Otherwise 
submit  three  copies  of  a  full-length,  camera-ready  paper  to  the  Technical  Program  Chairman.  Please  supply 
the  following  data  for  the  corresponding  authors:  name,  address,  email  address,  FAX,  and  phone  numbers. 

Authors  notified  of  acceptance  by  December  1, 1999. 


Microwave  components  Wireless  Radar  Imaging 

EMPEMI/EMC  Shielding  Fiberoptics  Radar  cross  section 

MIMIC  technology  Visualization  Fiberoptics 

Wave  propagation  Eddy  currents  Direction  finding 

Bioelectromagnetics  Inverse  scattering  Antennas 


PAPER  FORMATTING  REQUIREMENTS 

The  recommended  paper  length  is  6  pages,  with  8  pages  as  a  maximum,  including  figures.  The  paper  should 
be  camera-ready  (good  resolution,  clearly  readable  when  reduced  to  the  final  pnnt  of  6x9  mch  paper).  The 
paper  should  be  printed  on  8-l/2xl  1  inch  papers  with  13/16  side  margins,  1-1/16  inch  top  margin,  and  1  inch 
on  the  bottom.  On  the  first  page,  place  title  1-1/2  inches  from  top  with  authors,  affiliations  and  e-mail 
addresses  beneath  the  title.  Single  spaced  type  using  10  or  12  point  font  size,  entire  text  should  be  justified 
(flush  left  and  flush  right).  No  typed  page  numbers,  but  number  your  pages  lightly  in  pencil  on  the  back  ot 

each  page. 

For  all  questions  regarding  the  ACES  Symposium  please  contact: 

Douglas  H.  Werner,  Technical  Program  Chair 

The  Pennsylvania  State  University,  21 1A  Electrical  Engineering  East,  University  Park,  PA  16802 
Tel:  (814)  863-2946,  Fax:  (814)  865-7065,  E-mail:  aces@engr.psu.edu 
or  visit  ACES  on  line  at:  http://aces.ee.olemiss.edu  and  www.emclab.umr.edu/aces. 

EARLY  REGISTRATION  FEES 

ACES  member  $300  Student/Retired/Unemployed  $130  (no  proceedings) 

Non-member  $350  Student/Retired/Unemployed  $165  (includes  proceedings) 

Each  conference  registration  is  entitled  to  publish  two  papers  in  the  proceedings  free  of  charge.  Excess  pages 
over  a  paper  limit  of  8  will  be  charged  $  15/page. 

$500  BEST-PAPER  PRIZE 

A  $500  prize  will  be  awarded  to  the  authors  of  the  best  non-student  paper  accepted  for  the  16th  Annual 
Review.  Papers  will  be  judged  by  a  special  ACES  prize-paper  Committee  according  to  the  following 

1  Based  on  established  electromagnetic  (EM)  theory  4.  Practical  applications 

2  Reliable  data  5-  Estimates  of  computational  errors 

3!  Computational  EM  results  6.  Significant  new  conclusions 

$200  BEST  STUDENT  PAPER  CONTEST 

This  will  be  for  the  best  student  paper  accepted  for  the  16th  Annual  Review.  (Student  must  be  the  presenter 
on  the  paper  chosen).  Submissioiis  will  be  judged  by  three  (3)  members  of  the  BoD.  The  prizes  for  the 
student  presenter  and  his/her  principal  advisor  will  consist  of:  (1)  free  Annual  Review  regisfration  for  the 
following  year;  (2)  one  free  short  course  taken  during  the  2000  or  2001  Annual  Review;  and  (3)  $200  cash 

for  the  paper. 
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IEEE  Antennas  and  Propagation  Society  Russian  Foundation  for  Basis  Researches 


STEERING  COMMITTEE 

V.M.  Babich,  Steklov  Math.  Inst.,  St.  Petersburg 
V.S.  Buldyrev,  St.  Petersburg  Univ.,  St.  Petersburg 
G.  Pelosi,  Univ.  of  Florence,  Florence,  Italy 
J.L.  Volakis,  Univ.  of  Michigan,  Ann  Arbor,  USA 
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V.E.  Grikurov  (secretary) 

A.P.  Kiselev,  M.A.  Lyalinov 
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NEW  DIRECTIONS  IN  ASYMPTOTIC  TECHNIQUES 
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DIFFRACTION  ON  NON-SMOOTH  OBSTACLES 

EM  ANALYSIS  TECHNIQUES 
UNDERWATER  ACOUSTICS 
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Proposals  for  session’s  organizing  concerning  various  nature  wave  phenomena  are  encouraged! 


Working  language  is  English 


The  registration  fee  is  150  USD 


IMPORTANT  DEADLINES 

Proposals  for  session  organizing 
Papers  submissions 
Acceptance  of  papers 
Registration 


ASAP 

Feb.  1,2000 
March  1,  2000 
April  1,2000 


Early  registration  is  necessary  to  finalize  visa  formalities  in  time. 
Registration  fee  is  encouraged  but  not  requested  at  this  stage. 


WAYS  OF  SUBMISSION 


1.  Electronic  submission  (mostly  preferred!) 

I.V.  Andronov  iva@aa2628.spb.edu 

V.E.  Grikurov  orikurov@mph.phvs.spbu.ru 

A.P.  Kiselev  kiselev@pdmi.ras.ru 

2.  FAX:  +7-812-428-7240 

3.  Hard  copies  to  be  addressed:  Dept.  Math.&  Comp.Phys. 
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Ul’yanovskaya  -j  198904  Petrodvoretz 
Russia 

We  expect  to  receive  camera-ready  abstract  limited,  if  possible,  by  one  A4  page  of  lOpt  font  size 

TeX  and  Word  documents  are  accepted 

PUBLICATIONS 


Booklet  of  abstracts  of  all  papers  accepted  for  presentation 
will  be  available  at  the  Seminar. 

Full  manuscripts  will  be  invited  for  the  peer  selection 
to  be  included  into  either  the  Proceedings 
or  the  special  issue  of  IEEE  Antennas  and  Propagation  Transactions. 

PREREGISTRATION 

Pre-registration  is  necessary  to  be  included  into  mailing  list.  Please  fill  this  form  (all  fields  are 
mandatory)  and  return  it  to  one  of  addresses  given  above  as  soon  as  possible.  Updated 
information  and  on-line  pre-registration  available  via  URL  http://mph.phys.spbu.ru/DDMW 
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5TH  INTERNATIONAL  WORKSHOP  ON 


FINITE  ELEMENTS  FOR  MICROWAVE  ENGINEERING 
FINAL  CALL  FOR  PAPERS 
GRAND  CHALLENGES  FOR  NEW  MILLENNIUM 

John  Hancock  Conference  Center,  Boston,  Massachusetts,  USA 

June  8-9, 2000 

The  5th  Finite  Elements  Workshop  for  Microwave  Engineering  will  be  organized  by  Electrical  and  Computer  Engineering 
Department,  Worcester  Polytechnic  Institute  in  cooperation  with  University  of  Florence,  Italy  and  will  be  held  on  June 
8-9, 2000  at  John  Hancock  Conference  Center,  Boston,  Massachusetts,  U.S.A. 

The  workshop  provides  an  international  forum  for  reporting  and  discussing  recent  progresses  and  advances  in  the  finite 
elementtechnologies  for  microwave  engineering.  The  details  of  the  workshop  can  be  found  in  http://ece.wpi.edu~jinlee. 

February  1, 2000:  One-page  Abstract  due 

March  1,  2000:  Acceptance  Notification  will  be  mailed  to  the  corresponding  author  of  the  paper. 

April  15, 2000:  Pre-Register 

Authors  are  invited  to  submit  an  original  (camera-ready)  one-page  abstract  of  no  less  than  250  words.  The  abstract 
should  explain  clearly  the  content  and  relevance  of  the  proposed  contribution.  No  acknowledgments  should  be  included. 
Your  cover  letter  should  includethe  complete  mailing  address,  telephone,  fax  number,  and  e-mail  address  (if  available) 
for  the  corresponding  author.  Please  mail  abstracts-do  not  send  via  facsimile.  Fulll  paper  is  due  at  the  conference  and 
will  be  under  regular  peer  review  process.  Accepted  papers  will  be  published  in  a  special  issue  of  Electromagnetics, 
2001. 

Each  presenting  author  is  required  to  register,  via  a  non-refundable  pre-registration  fee  of  US$1 00  and  must  be  sent 
by  April  15,  2000.  The  registration  fee,  to  be  determined,  will  include  the  Workshop  program  and  proceedings, 
attendance  at  all  technical  sessions,  refreshments,  the  opening  reception,  and  the  banquet. 

Local  arrangements  should  be  directed  to:  Abstracts  should  be  directed  to: 

Ms.  Yurong  Sun  Prof.  Jin-Fa  Lee 

5th  Finite  Element  Wkshop  for  Microwave  Engineering  5th  Finite  Element  Wkshop  for  Microwave  Engineering 
EM/CAD  Lab,  ECE  Department,  WPI  EMCAD  Lab,  ECE  Dept,  WPI. 

1 00  Institute  Road  1 00  Institute  Road 

Worcester,  MA  01609,  USA  Worcester,  MA01609,  USA 

Tel:  508-831-5757  Fax:  508-831-5491  Tel:  508-831-5778,  Fax:  508-831-5491 

Email:  ysun@ece.wpi.edu  Email:  jinlee@ece.wpi.edu 

ABSTRACT  GUIDELINES 

Submit  in  English.  Use  single-spaced,  8.5x1 1  inch  paper;  using  1 2  point  Times  Roman  or  and  equivalent  serif  typeface. 
Set  margins  to  25mm  (1  inch)  and  set  paragraph  indentation  to  3.5mm  (0.14  inch).  Type  title  in  bold  letters;  centered 
at  the  top.  Below  title,  center  name  of  the  author(s)  with  affiliation  and  complete  mailing  address. 

ADVANCE  PROGRAM,  AND  INFORMATION  SENT  BY  MAY  1 , 2000 

Registration:  Travel,  lodging,  registration  and  local  information  will  be  mailed  with  Advance  Program  by  May  1 , 2000. 
The  advance  registration  fee  for  all  participants,  including  session  chairpersons  and  authors,  is  US$100.  The  fee 
includes  the  Workshop  program  and  proceedings,  attendance  at  all  technical  sessions,  refreshments,  the  opening 
reception  and  the  banquet. 


LOCATION 


The  Workshop  is  planned  to  be  held  at  the  John  Hancock  Conference  Center,  Boston,  MA,  USA,  on  Thursday  and 
Friday,  June  8-9,  2000.  Boston  is  the  favorite  conference  venue  in  the  Northeast  US  and  one  of  the  most  popular 
conference  cities  in  the  US.  The  city  is  beautiful  and  alive.  The  charm  and  attractions  of  historic  Boston  should  create 
an  enjoyable  environ  ment  for  the  workshop. 
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limitation  removal,  or  other  performance  improvement. 
Note:  Code  (or  algorithm)  capability  descriptions  are 
not  acceptable,  unless  they  contain  sufficient  technical 
material  to  justify  consideration. 

7.  Code  input/output  issues.  This  normally  involves 
innovations  in  input  (such  as  input  geometry 
standardization,  automatic  mesh  generation,  or  computer- 
aided  design)  or  in  output  (whether  it  be  tabular,  graphical, 
statistical,  Fourier-transformed,  or  otherwise  signal- 
processed).  Material  dealing  with  input/output  database 
management,  output  interpretation,  or  other  input/output 
issues  will  also  be  considered  for  publication. 

8.  Computer  hardware  issues.  This  is  the  category  for 
analysis  of  hardware  capabilities  and  limitations  in  meeting 
various  types  of  electromagnetics  computational  require¬ 
ments.  Vector  and  parallel  computational  techniques  and 
implementation  are  of  particular  interest. 


Applications  of  interest  include,  but  are  not  limited  to, 
antennas  (and  their  electromagnetic  environments), 
networks,  static  fields,  radar  cross  section,  shielding, 
radiation  hazards,  biological  effects,  electromagnetic  pulse 
(EMP),  electromagnetic  interference  (EMI),  electromagnet¬ 
ic  compatibility  (EMC),  power  transmission,  charge 
transport,  dielectric  and  magnetic  materials,  microwave 
components,  MMIC  technology,  remote  sensing  and  geo¬ 
physics,  communications  systems,  fiber  optics,  plasmas, 
particle  accelerators,  generators  and  motors,  electromagnet¬ 
ic  wave  propagation,  non-destructive  evaluation,  eddy 
currents,  and  inverse  scattering. 

Techniques  of  interest  include  frequency- domain  and 
time-domain  techniques,  integral  equation  and  differential 
equation  techniques,  diffraction  theories,  physical  optics, 
moment  methods,  finite  differences  and  finite  element 
techniques,  modal  expansions,  perturbation  methods,  and 
hybrid  methods.  This  list  is  not  exhaustive. 

A  unique  feature  of  the  Journal  is  the  publication  of 
unsuccessful  efforts  in  applied  computational 
electromagnetics.  Publication  of  such  material  provides  a 
means  to  discuss  problem  areas  in  electromagnetic  model¬ 
ing.  Material  representing  an  unsuccessful  application  or 
negative  results  in  computational  electromagnetics  will  be 
considered  for  publication  only  if  a  reasonable  expectation 
of  success  (and  a  reasonable  effort)  are  reflected. 
Moreover,  such  material  must  represent  a  problem  area  of 
potential  interest  to  the  ACES  membership. 

Where  possible  and  appropriate,  authors  are  required  to 
provide  statements  of  quantitative  accuracy  for  measured 
and/or  computed  data.  This  issue  is  discussed  in  "Accuracy 
&  Publication:  Requiring  quantitative  accuracy  statements 
to  accompany  data",  by  E.K.  Miller,  ACES  Newsletter ,  Vol. 
9,  No.  3,  pp.  23-29,  1994,  ISBN  1056-9170. 

EDITORIAL  REVIEW 

In  order  to  ensure  an  appropriate  level  of  quality  control, 
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the  computer  codes,  computational  techniques,  and 
applications  discussed  in  the  paper  (as  applicable)  and 
should  otherwise  be  usable  by  technical  abstracting  and 
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3.  Either  British  English  or  American  English  spellings 
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consistently  throughout  the  paper. 

4.  Any  commonly-accepted  format  for  referencing  is 
permitted,  provided  that  internal  consistency  of  format  is 
maintained.  As  a  guideline  for  authors  who  have  no  other 
preference,  we  recommend  that  references  be  given  by 
author(s)  name  and  year  in  the  body  of  the  paper  (with 
alphabetical  listing  of  all  references  at  the  end  of  the  paper). 
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Titles  of  papers  or  articles  should  be  in  quotation  marks. 

5.  Internal  consistency  shall  also  be  maintained  for  other 
elements  of  style,  such  as  equation  numbering.  As  a 
guideline  for  authors  who  have  no  other  preference,  we 
suggest  that  equation  numbers  be  placed  in  parentheses  at 
the  right  column  margin. 

6.  The  intent  and  meaning  of  all  text  must  be  clear.  For 
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material  is  neat,  legible,  and  reproducible.  The  preferred 
font  style  is  Times  Roman  10  point  (or  equivalent)  such  as 
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that  used  here  is  preferred.  No  author's  work  will  be 
turned  down  once  it  has  been  accepted  because  of  an 
inability  to  meet  the  requirements  concerning  fonts  and 
format.  Full  details  are  sent  to  the  author(s)  with  the  letter 
of  acceptance. 
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over  "nth-generation"  photocopies.  These  original  figures 
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